DEVELOPMENT  OF  A 
HIGHER-ORDER  UPWIND  ALGORITHM 
FOR  DISCONTINUOUS  COMPRESSIBLE  FLOW 

THESIS 


Barry  A.  Croker,  1LT,  USAF 
AFIT/G  AE /ENY/05-M05 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the  United 
States  Government. 


AFIT/GAE/ENY/05-M05 


DEVELOPMENT  OF  A 
HIGHER-ORDER  UPWIND  ALGORITHM 
FOR  DISCONTINUOUS  COMPRESSIBLE  FLOW 


THESIS 


Presented  to  the  Faculty 
Department  of  Aeronautics  and  Astronautics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Aeronautical  Engineering 


Barry  A.  Croker,  B.S.M.E 
1LT,  USAF 


March  2005 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AFIT/GAE/ENY/05-M05 


DEVELOPMENT  OF  A 
HIGHER-ORDER  UPWIND  ALGORITHM 
FOR  DISCONTINUOUS  COMPRESSIBLE  FLOW 


Barry  A.  Croker,  B.S.M.E 
1LT,  USAF 


Approved: 


/signed/ 


21  Mar  2005 


Maj  Richard  J.  McMullan,  Ph.D.  date 

(Chairman) 


/signed/ 


21  Mar  2005 


Lt  Col  Raymond  C.  Maple,  Ph.D.  date 

(Member) 


/signed/ 


21  Mar  2005 


Lt  Col  Eric  J.  Stephen,  Ph.D. 
(Member) 


date 


AFIT/GAE/ENY/05-M05 


Abstract 

The  use  of  compact  differences  to  discretize  partial  differential  equations  produces 
higher  orders  of  accuracy  with  a  smaller  stencil  and  lower  absolute  error  than  approxima¬ 
tions  using  traditional  finite  differences.  However,  most  solution  algorithms  for  the  com¬ 
pressible  Euler  equations  that  include  compact  differencing  cannot  be  applied  to  problems 
involving  discontinuities  (such  as  shock  or  contact  surfaces)  without  the  use  of  filtering  or 
limiting  methods.  A  global  fourth-order  method  that  incorporates  compact  differencing  with 
Roe’s  approximate  Riemann  solver  was  investigated.  This  method  was  incorporated  into  a 
one-dinrensional  numerical  simulation  of  the  compressible  Euler  equations,  and  applied  to  a 
one-dinrensional  shock  tube  problem.  The  method  was  also  extended  to  two  dimensions,  and 
applied  to  a  two-dimensional  shock  tube  problem  and  an  advecting  vortical  structure  prob¬ 
lem  on  both  rectilinear  and  curvilinear  meshes.  The  results  were  compared  to  a  third-order 
Roe  scheme  and  a  fourth-order  compact  difference  scheme.  The  compact  Roe  scheme  was 
found  to  be  consistent  in  both  the  one-dinrensional  and  two-dimensional  simulations.  An 
order  of  accuracy  determination  showed  that  it  has  an  order  of  accuracy  somewhere  near 
fourth  order,  with  absolute  error  comparable  to  that  of  the  standard  compact  difference 
scheme.  With  proper  selection  of  solution  parameters,  the  scheme  was  also  shown  to  ac¬ 
curately  capture  a  discontinuous  solution  where  unfiltered  compact  schemes  would  become 
unstable. 
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DEVELOPMENT  OF  A 


HIGHER-ORDER  UPWIND  ALGORITHM 
FOR  DISCONTINUOUS  COMPRESSIBLE  FLOW 

I.  Introduction 

1.1  Modeling  The  Governing  Equations  of  Fluid  Flow 

The  governing  equations  of  fluid  flow  are  hyperbolic,  non-linear,  partial  differential 
equations  (PDE)  of  multiple  variables.  They  are  continuum  statements  of  conservation  of 
mass,  momentum  and  energy,  which  can  describe  the  flow  of  any  Newtonian  fluid  under 
almost  any  circumstance.  They  were  originally  proposed  by  Leonhard  Euler  in  the  mid 
1750’s,  and  later  refined  by  Claude  Louis  Marie  Henri  Navier  in  the  early  1820’s.  Nearly 
two  hundred  years  later,  their  work  still  remains  somewhat  unfinished.  Because  of  the 
complexity  of  these  governing  equations,  no  generalized  closed-form,  exact  solutions  have 
been  discovered  [I]. 

For  years  scientists  and  engineers  relied  on  various  methods  to  approximate  the  so¬ 
lutions  to  these  equations.  Many  of  these  solutions  were  limited  to  specialized  or  idealized 
cases,  with  simplifications  necessary  to  bring  the  problem  closure.  The  advent  of  the  com¬ 
puter  age  opened  a  new  chapter  in  the  field  of  solving  these  differential  equations.  No  longer 
did  solutions  have  to  be  simplified  and  terms  ignored  to  be  solved.  Computers  allowed  the 
governing  equations  to  be  approximated  through  various  linearization  and  discritization 
techniques  and  to  be  solved  numerically. 

As  one  can  imagine,  these  numerical  solution  techniques  do  not  involve  trivial  compu¬ 
tations.  Current  computer  technology  has  yet  to  advance  to  the  point  that  many  complex 
flow  problems  can  be  fully  modeled  and  accurately  solved.  However,  the  technology  to  solve 
what  was  once  a  daunting  problem  for  a  fluid  engineer  can  now  be  approximated  quite 
accurately  by  anyone  who  possesses  a  personal  computer. 

In  their  most  general  form,  the  governing  equations  have  two  major  components — 
partial  derivatives  with  respect  to  time  (temporal  derivatives),  and  partial  derivatives  with 
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respect  to  space  (spatial  derivatives).  The  temporal  derivatives  are  somewhat  straightfor¬ 
ward  to  approximate,  as  time  only  increases  in  one  dimension  and  in  an  orderly  fashion.  The 
real  challenge  is  in  developing  an  accurate  approximation  to  the  spatial  derivative  terms. 
The  methods  chosen  to  approximate  the  spatial  derivatives  are  highly  dependent  on  the 
problem  being  solved.  They  are  influenced  by  the  problem  geometry,  the  flow  conditions, 
and  even  the  fluids  involved. 


1.2  Discretization  of  the  Governing  Equations 

These  governing  equations  are  continuum  statements,  which  means  they  can  describe 
physical  phenomena  at  any  point  in  space  at  any  time.  The  first  step  in  trying  to  approx¬ 
imate  them  is  to  discretize  them,  or  divide  them  up  into  finite  portions  in  space  and  time. 
The  value  of  a  continuous  function  at  a  discrete  point  in  space  can  be  approximated  by 
the  value  of  the  function  and  its  derivatives  at  other  discrete  points.  This  approximation  is 
called  a  Taylor  Series  Expansion  (TSE).  For  example,  the  TSE  of  a  function  /  at  (x0  + Ax) 
can  be  expressed  as  a  function  of  the  value  at  xQ  and  the  derivatives  of  /  at  x0. 


f{xQ  +  Ax)  =  f(xQ)  + 


df 

dx 


Ax  + 


S2/ 

dx2 


Ax2  d3f 
+  dx 3 


Ax3  dnf 

- (-...-( - - 

3!  dxn 


Ax’1 


ni 


(1.1) 


Because  of  the  factorial  terms  in  the  denominator,  a  TSE  is  a  converging  series.  If 
the  TSE  included  an  infinite  number  of  terms  it  would  produce  the  exact  value  of  /.  Since 
only  a  finite  number  of  terms  can  actually  be  used,  the  approximation  has  an  error  term 
associated  with  it.  The  magnitude  of  the  exponent  on  Ax  in  the  last  term  of  the  TSE  is 
called  the  order  of  the  expansion. 


1.3  Finite  Difference  and  Compact  Difference  Approximations 

If  the  terms  in  the  TSE  are  re-arranged,  discrete  approximations  for  the  derivative 
terms  can  be  developed.  The  exact  value  of  the  derivative  can  be  expressed  as  the  sum 
of  a  finite-difference  expression  and  a  truncation  error.  For  example,  if  the  first  derivative 
term  is  moved  to  the  left-hand  side  and  solved  for,  the  result  is  an  expression  of  the  first 
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derivative  df/dx  as  a  function  of  the  values  of  /  at  two  known  points: 


df_  =  f(x0  +  Ax)  -  f(x0)  8Pj_  Ax 
dx  Ax  dx2  2! 


The  finite-difference  expression  is  the  first  term  on  the  left-hand  side,  while  the  truncation 
error  is  an  infinite  series  beginning  with  the  d2  f  /dx2  term.  This  finite-difference  approxi¬ 
mation  is  first-order,  as  Ax  in  the  leading  truncation  error  term  has  an  exponent  of  one.  If 
terms  in  the  truncation  error  are  replaced  with  other  finite-difference  expressions,  higher  ac¬ 
curacy  approximations  can  be  derived.  For  example,  a  fourth  order  accurate  representation 
of  the  first  derivative  of  the  variable  /  at  point  i  could  be  expressed  as: 

dfi  _  fj- 2  -  8/i-i  +  8/i+i  -  fj+2  ^  A  ,4 
dx  12Ax  X 


In  essence,  the  derivative  at  a  point  is  determined  by  a  curve  fit  through  the  values  at 
surrounding  points.  In  this  case,  the  finite  difference  expression  is  a  central  difference  where 
information  from  either  side  is  used  to  determine  the  slope  at  the  point  in  question.  Finite 
differences  can  also  be  one  sided,  either  forwards  or  backwards,  using  information  from  only 
one  direction. 

Since  the  governing  PDEs  consist  of  multiple  partial  derivatives,  these  discretized 
derivatives  can  be  substituted  and  numerical  approximations  made.  Additional  work  is 
needed  for  the  linearization  of  certain  terms,  but  the  underlying  methodology  remains  the 
same:  continuous  functions  are  replaced  with  discretized  representations,  each  having  their 
associated  truncation  error.  The  combined  influence  of  the  individual  truncation  errors 
results  in  an  overall  order  of  accuracy  of  the  solution. 

Higher  order  solutions  result  in  better  approximations  due  to  both  the  magnitude  of 
the  truncation  error  and  the  rate  at  which  the  truncation  error  approaches  zero.  Since  the 
truncation  error  is  a  function  of  Ax,  the  accuracy  of  the  solution  is  often  dependent  on 
the  refinement  of  the  mesh.  Since  most  solution  meshes  have  values  of  Ax  <  1,  the  higher 
the  order  the  approximation  the  smaller  the  truncation  error  terms.  When  comparing 
two  solutions  of  different  orders  a  higher  order  approximation  allows  for  better  resolution 
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and  accuracy  of  flow  structures  and  behaviors  without  the  need  of  a  greatly  refined  grid. 
Conversely,  a  low  order  approximation  on  a  coarse  grid  will  be  more  prone  to  error  than  a 
higher  order  one  would  be  on  the  same  grid. 

In  recent  years,  advances  in  computing  power  have  allowed  methods  such  as  large 
eddy  simulation  (LES)  and  direct  numerical  simulation  (DNS)  to  gain  popularity.  Both  of 
these  methods  attempt  to  numerically  solve  the  governing  equations  of  turbulent  fluid  flow 
completely  from  large  to  very  small  scales.  Since  many  turbulent  problems  involve  complex 
flow  patterns  and  very  small  length  scales,  these  methods  require  considerable  computing 
power  and  highly  refined  meshes.  First  or  second  order  approximations  often  cannot  resolve 
the  high  frequencies  and  miniscule  length  scales  involved  in  these  types  of  simulations. 
Information  propagating  at  the  high  frequencies  is  overcome  by  the  truncation  error  terms, 
and  the  solution  begins  to  break  down.  Small  scale  flow  structures  are  dissipated  out  and 
the  turbulent  structures  are  lost. 

Traditional  finite  differences  can  in  theory  be  derived  for  any  accuracy,  but  in  practice 
they  are  limited  to  second  or  third  order.  The  order  of  a  finite  difference  is  related  to  the 
number  of  points  included  in  the  “stencil”,  or  function  points  used  in  the  calculation.  This 
number  of  points  also  increases  with  the  degree  of  the  derivative.  For  example,  for  a  given 
accuracy  a  second  derivative  needs  more  points  than  a  first  derivative.  A  wider  stencil 
not  only  requires  a  more  complicated  solution  algorithm  but  can  cause  complications  at 
boundaries  of  the  solution  domain. 

An  alternate  approach  to  gaining  higher  order  approximations  without  as  large  a  sten¬ 
cil  size  is  through  the  use  of  Pade,  or  compact,  finite  difference  methods  jl5].  In  the  compact 
formulation,  the  value  of  the  derivative  f  is  still  expressed  as  a  function  of  the  surrounding 
values  of  /,  but  in  an  implicit  relation.  For  example,  a  fourth-order  representation  of  the 
first  derivative  f  can  be  expressed  as 

\fU  +  fi  +  \fU  1  =  \  <L4> 

In  this  case  the  unknowns  are  the  f  terms  on  the  left-hand  side  of  the  equation,  while  the 
known  values  are  on  the  right-hand  side.  For  a  domain  1  <  i  <  n,  a  system  of  n  equations 
can  be  written  and  solved  simultaneously  for  n  unknown  values  of  f .  In  this  case  the 
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number  of  points  in  the  approximation  has  been  reduced  from  five  points  to  three,  making 
the  stencil  more  “compact” . 

1-4  Advantages  of  the  Compact  Difference 

One  major  advantage  of  using  compact  difference  methods,  aside  from  having  a  smaller 
stencil,  is  an  increase  in  solution  resolution  over  traditional  finite  difference  methods.  Infor¬ 
mation  in  a  hyperbolic  domain  propagates  via  waves  of  characteristics.  A  Fourier  decompo¬ 
sition  of  a  solution  would  show  a  broad  spectrum  of  waves  from  large  to  small  wavenumbers. 
Certain  portions  of  the  flow  structure  rely  on  different  portions  of  this  spectrum  to  propagate 
information.  Large  scale  eddys  transport  energy  through  the  low  frequency  waves,  while 
turbulent  scales  dissipate  energy  at  the  high  frequencies.  As  previously  stated,  it  is  these 
high-frequency  waves  that  are  the  most  challenging  to  model.  Lele  (T5J  defines  a  resolving 
efficiency  as  the  ratio  of  the  maximum  wave  number  for  a  well  resolved  wave  (|/cAx|  <  n) 
to  the  theoretically  existing  Fourier  wavenumbers.  Whereas  traditional  fourth-order  differ¬ 
ences  have  a  resolving  efficiency  near  e  =  0.19,  a  compact  formulation  efficiency  is  closer  to 
e  =  0.477  (l2  .  This  means  nearly  twice  as  much  high-frequency  information  is  captured  by 
the  compact- difference  as  by  the  finite-difference.  Additionally,  compact  difference  meth¬ 
ods  also  exhibit  smaller  phase  speed  errors  further  resulting  in  better  resolution  of  high 
wavenumber  flow  features.  Both  of  these  features  mean  that  an  LES  or  DNS  solution  using 
compact-differencing  will  almost  always  produce  a  more  accurate  solution  than  one  using 
traditional  finite-differences. 

Another  consequence  of  using  compact  difference  methods  is  a  decrease  in  the  mag¬ 
nitude  of  the  truncation  error,  and  hence  an  additional  increase  in  the  overall  accuracy  of 
the  solution.  While  a  finite-difference  and  a  compact- difference  scheme  could  both  have  a 
truncation  error  term  of  the  same  order,  the  magnitude  of  this  term  is  generally  smaller 
for  the  compact  case.  Because  of  this  difference  in  absolute  error,  a  fourth  order  compact 
scheme  of  the  first  derivative  can  have  a  Ax  of  63  «  1.57  larger  than  a  standard  fourth 
order  finite  difference  expression  while  maintaining  the  same  accuracy  |l6|. 
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1.5  Problems  with  Compact  Differencing 


With  better  high  wavenumber  resolution,  lower  overall  error,  and  a  smaller  stencil,  the 
compact  difference  at  first  appears  to  be  a  superior  method  for  spatial  discretization.  For  the 
approximation  of  a  linear  function  this  is  usually  the  case,  but  one  must  remember  the  PDE’s 
which  govern  generalized  fluid  flow  are  highly  non-linear,  and  additional  considerations 

uses 


must  be  taken.  The  right-hand  side  of  a  compact  difference  such  as  equation  1.4 


information  from  before  and  after  the  desired  evaluation  point,  i.e.  for  f[  the  values  of 
fi-i  and  /*_)_ i  are  needed.  This  “before  and  after”  is  termed  a  central  difference,  in  this 
case  centered  at  i.  For  non-linear  problems  which  contain  discontinuous  solutions,  central 
difference  methods  are  unconditionally  unstable  and  will  oscillate  without  bound.  This  is 
caused  by  the  fact  that  the  characteristic  information  only  propagates  in  a  single  direction. 
If  information  is  incorporated  from  a  location  that  is  downstream  of  the  direction  of  the 
characterics,  non-physical  behavior  will  result.  Since  the  compact-difference  scheme  uses 
a  central  difference  term  on  the  right-hand  side,  it  is  unusable  for  problems  containing 
discontinuities,  specifically  problems  involving  shock  waves. 


One  method  which  allows  compact  differences  to  be  used  with  such  problems  is  to 
use  spectral- like  filtering  methods  to  remove  non-physical  oscillations.  These  spectral  filters 
remove  the  unbounded  high-frequency  oscillations  before  they  have  a  negative  impact  on  the 
overall  solution.  To  be  efficient,  a  spectral  filter  should  have  twice  the  order  of  accuracy  as 
that  of  the  solution  method  ||5j.  For  example,  on  a  fourth-order  solution,  an  eight  or  higher 
order  filter  would  be  required.  Not  only  do  these  high-order  filters  require  a  large  stencil, 
they  are  computationally  expensive  and  add  considerably  to  the  complexity  of  the  solution. 
An  alternate  method,  which  is  the  goal  of  this  thesis,  is  to  incorporate  an  upwind  bias  into 
the  compact  difference  formulation  and  eliminate  the  sources  of  oscillations  altogether. 

An  upwind  bias  means  that  the  solution  method  only  uses  information  that  is  relevant — 
information  that  is  “upwind”  of  the  point  in  question.  An  upwind  bias  for  the  previous 
example  could  be  the  value  of  f-  would  depend  on  /*_ i  and  /*_ 2,  while  not  considering  fj+  \ . 
Even  with  the  use  of  standard  finite  differences,  upwind  biasing  is  necessary  for  an  accurate, 
physical  solution.  If  upwinding  could  be  incorporated  into  a  compact  difference  scheme,  a 
global  higher-order  solution  could  be  developed  that  would  forgo  the  use  of  filters.  Work 
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by  Ravichandran  |T8]  demonstrated  that  it  is  feasible  to  incorporate  the  upwinding  into 
the  compact  formulation  using  flux  vector  splitting,  where  the  sign  of  the  flux  Jacobians 
determines  the  direction  of  information  propagation.  However,  the  derived  scheme  in  itself 
is  not  non-oscillatory,  and  requires  the  use  of  limiters  at  solution  extrema. 


1.6  Development  of  an  Upwind  Compact- Difference  Formulation 

It  was  suggested  by  David  Gottlieb  [7]  that  using  an  upwind  flux- differencing  scheme, 


such  as  Roe’s  scheme,  directly  into  a  compact  formulation  such  as  Equation  1.4  could 


allow  for  a  stable  higher  order  solution  method  without  the  use  of  filters.  Fluxes  would 
be  calculated  using  Roe’s  scheme,  then  substituted  directly  into  the  compact  scheme  to 
produce  a  higher  order  approximation  of  df/dx.  Since  many  existing  numerical  simulations 
include  a  Roe  flux  solver,  this  could  imply  that  higher-order  solutions  could  be  calculated 
readily  without  a  new  solution  algorithm.  This  proposition  could  be  summarized  as  follows: 


aH-i  +  F'i  +  aH+i  ~  ^(-^i+i/2  ~  Fi- 1/2) 
where  F  denotes  Roe  averaged  fluxes,  and 


(1.5) 


«=g(3-2«) 


a  =; 


22 

Preliminary  investigation  shows  that  this  solution  formulation  is  only  partially  correct. 


While  the  once  unconditionally- unstable  compact  scheme  of  Equation  1.4  is  now  made 
stable,  it  appears  to  only  have  overall  order  the  same  as  the  order  of  the  flux  calculations. 
Hence,  if  F  is  calculated  to  first  order  the  compact  formulation  will  only  produce  first 
order  F'  values.  Deng  and  Maekawa  |3|  have  developed  a  fourth-order  compact  variable 


interpolation  scheme  which  appears  to  allow  method  Equation  1.5  to  be  globally  fourth 
order.  Their  findings  only  present  results  in  one-dimension  and  attempt  no  verification  of 
the  formal  order  of  accuracy  of  the  new  method. 

It  is  the  goal  of  this  thesis  to  investigate  the  proposition  of  Gottlieb  and  to  better 
understand  the  new  method  of  Deng  and  Maekawa  |3|.  This  new  “cRoe”  method  is  based 
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on  compact- difference  methods  and  Roe’s  approximate  Riemann  solver,  but  forgoes  the 
use  of  spectral  filters  through  the  use  of  an  adaptive  compact  variable  interpolant.  The 
development  of  the  regular  compact  difference  will  first  be  discussed.  This  compact  method 
will  then  be  applied  to  a  model  linear  equation  in  one-dinrension  and  compared  with  a 
first-order  explicit  method.  Next,  the  cRoe  method  will  be  applied  to  a  one-dinrensional 
discontinuous  Riemann  problem  and  compared  with  a  third-order  Roe  scheme  (oRoe).  The 
cRoe  scheme  will  then  be  extended  to  two- dimensions  and  applied  to  a  similar  Riemann 
problem.  These  Riemann  problems  will  allow  for  a  measure  of  the  consistency  of  the  hybrid 
cRoe  scheme  while  also  allowing  a  comparison  with  the  oRoe  method.  Next,  the  cRoe 
method  will  be  compared  with  the  third-order  oRoe  and  a  fourth-order  compact  differencing 
scheme  (cC)  on  a  problem  involving  the  advection  of  a  vortical  structure  on  a  rectilinear 
mesh.  This  problem  will  allow  the  formal  order  of  accuracy  of  the  cRoe  solution  to  be 
determined  and  the  absolute  error  to  be  compared  to  the  oRoe  and  cC  schemes.  Finally, 
the  three  methods  will  be  applied  to  the  same  vortex  problem  on  a  curvilinear  mesh  to  gain 
a  more  realistic  measure  of  the  usefulness  of  the  solution  methods. 
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II.  Governing  Equations 

In  order  to  understand  the  behavior  of  various  schemes  involving  both  compact  and  finite 
differences  several  model  problems  will  be  considered.  To  verify  the  validity  of  the  compact 
difference  formulation  a  linear  one- dimensional  problem  will  first  be  investigated. 


2.1  ID  Linear  Wave  Equation 

The  1-D  linear  wave  equation  (LWE)  |2l]: 


df  ,  df_ 

Tt+aTx~° 


(2.1) 


(or  advection  equation)  is  a  hyperbolic  partial  differential  equation  that  describes  the 
motion  of  a  wave  moving  in  the  x  direction  at  a  velocity  a.  It  also  can  describe  the  convective 
transport  of  a  scalar  quantity  /  at  a  flow  speed  a  (TT] .  It’s  exact  solution  is  described  by 


f(x,  t )  =  F(x  —  aAt) 


(2.2) 


given  initial  data 


f(x,  0)  =  F(x)  —  oo  <  x  <  oo 


(2.3) 


The  LWE  is  a  useful  tool  which  can  model  the  transportive  behavior  of  the  equations  that 
govern  generalized  fluid  flow.  Since  the  LWE  is  linear  and  scalar  it  is  much  simpler  to 
use  than  fully  non-linear  vectorized  equations.  Solution  methods  for  nonlinear  equations 
are  often  initially  applied  to  the  LWE  as  a  preliminary  form  of  validation.  If  a  solution  is 
unstable  when  applied  to  the  LWE  it  will  almost  certainly  be  unstable  on  a  more  complex 
problem.  However,  if  a  solution  is  stable  when  applied  to  the  LWE  it  could  be  worth 
investigating  with  a  more  complex  equation  set  such  as  Eulers  equations. 


2.2  ID  Euler  Equations 

The  Euler  Equations  are  the  governing  equations  for  compressible  inviscid  fluid  flow. 
They  consist  of  two  scalar  equations  for  conservation  of  mass  and  energy  and  one  vector 
equation  for  conservation  of  momentum.  All  three  equations  form  a  first  order  system  of 
coupled  non-linear  partial  differential  equations  (PDEs). 
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The  Euler  equations  can  be  written  in  a  number  of  equivalent  forms  based  on  the 
derivation  method  used.  If  a  PDE  can  be  expressed  having  constant  coefficients  on  the 
derivative  terms,  (or  if  the  coefficients  are  variable  the  derivatives  appear  only  once),  the 
formulation  is  said  to  be  in  “conservation  form”.  “Divergence  form”  is  also  an  acceptable 
description  for  a  physical  conservation  law  with  divergence  terms  of  physical  quantities.  [2l]. 
In  order  to  compute  a  correct  solution  for  a  discontinuous  flow  the  conservation  form  of  the 
governing  equations  must  be  used  |ll] . 

Two  sets  of  variables  are  be  used  in  the  conservation  form  of  the  Euler  equations: 
conservative  variables  and  primitive  variables.  Conservative  variables  are  the  conserved 
quantities  that  the  equations  are  derived  from:  conservation  of  mass,  momentum  and  energy. 
Primitive  variables  are  those  that  can  be  experimentally  measured  and  are  often  imposed 
by  the  physical  domain  of  the  problem,  i.e.  density  p,  velocity  u,  and  pressure  p.  Written 
in  one-dinrensional,  conservative,  differential  form  the  Euler  set  is  given  as  [ll] : 


dU  dF  _ 
dt  dx 


(2.4) 


The  vector  of  conserved  variables  U  and  vector  of  conserved  fluxes  F  can  be  expressed  in 
terms  of  the  primitive  variables  as 


p 

pu 

u  = 

pu 

,F  = 

pu 2  +  p 

Et 

(- Et+p)u 

Equation  2.4  is  a  system  of  three  equations  with  four  unknowns:  p ,  u.  p,  and  Ef.  To  reach 
closure  an  additional  equation  of  state  is  added  which  relates  the  thermodynamic  properties 
of  the  fluid.  One  common  assumption  is  that  of  a  calorically  perfect  gas  such  that 


p  =  p1ZT 


(2.6) 
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where  7Z  is  the  universal  gas  constant  and  T  is  the  temperature.  With  this  assumption 
total  energy  can  be  calculated  using 


E,= 


P  ,  1  2 

- 1"  +  opu 

7  —  1  2 


(2.7) 


where  7  is  the  ratio  of  specific  heats.  For  this  thesis  the  fluid  considered  will  be  air  with 
7=1.4  and  71=287  N  rn/kg  K. 


The  Equation  2.4  can  be  used  directly  to  solve  the  Euler  set  using  a  finite  difference 
solution  method.  A  finite-difference  solution  method  computes  the  solution  to  a  flow  by 
looking  at  the  changes  in  the  conserved  variables  at  discrete  points  in  space  and  time.  The 
system  of  PDE’s  described  by  Euler’s  equations  represent  continuum  solutions  which  at  any 
point  in  time  or  space  represent  a  statement  of  conservation.  In  order  to  reach  a  solution 
to  these  equations  finite- difference  expressions  can  be  made  that  well  approximate  the  PDE 
at  a  localized  point.  If  these  approximations  individually  represent  the  actual  solution  it  is 
expected  that  the  global  solution  will  also  closely  match  the  exact  solution  and  the  method 
is  said  to  have  a  conservative  property  [2l].  As  more  and  more  points  are  added  to  the 
solution  method  these  approximations  become  closer  and  closer  to  the  exact  value  of  the 
PDE  and  solution  method  is  said  to  be  consistent  1 21  . 


An  alternate  solution  method  can  be  derived  from  looking  not  only  at  the  PDE  as 
a  mathematical  expression  but  as  a  description  of  a  physical  phenomena.  If  an  arbitrary 
control  volume  is  taken  in  the  solution  domain  it  can  be  shown  that  the  time  rate  of  change 
of  the  conserved  variables  within  the  control  volume  is  equal  to  the  net  flux  through  the 
control  volume  boundary  (Tlj.  In  this  method  the  Euler  equations  in  integral  form  are 


4-  [  UcN  +  [  Fi-dS  =  0 
dt  Jy  Js 


(2.8) 


where  U  and  F  have  the  same  definitions  as  the  previous  form.  This  result  could  also 


be  reached  by  taking  the  differential  form  of  Equation  2.4  and  applying  the  divergence 
theorrn  111  . 


A  finite-volume  solution  method  uses  this  integral  form  to  divide  the  computational 
domain  into  a  mesh  of  discrete  cells.  As  the  size  of  the  control  volumes  are  reduced,  i.e  the 
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mesh  is  refined,  it  is  expected  that  the  approximate  value  in  the  control  volume  approaches 
it’s  exact  value  and  the  method  is  again  considered  to  be  consistent. 


In  order  for  Equation  2.8  to  be  conservative,  it  must  hold  true  for  any  sized  control 


volume.  This  means  that  a  finite-volume  solution  method  does  not  require  a  uniform  com¬ 
putational  domain.  Since  most  physical  flow  problems  cannot  be  easily  discretized  with  a 
uniform  mesh  a  finite  volume  method  is  useful  for  “real  world”  problems  1 11  . 


2.3  2D  Euler  Equations 

When  extended  to  two  dimensions  the  Euler  set  changes  slightly  due  to  an  additional 
degree  of  freedom.  Whereas  in  one  dimension  there  was  only  a  single  velocity,  in  two 
dimensions  there  are  two  separate  velocity  components:  u  in  the  x  direction  and  v  in  the 
y  direction.  Combined  with  mass  these  two  velocities  make  up  two  separate  momentums 
which  must  be  individually  conserved.  Additionally,  whenever  a  velocity  is  considered  in  an 
energy  term  the  total  vector  sum  of  u  and  v  must  now  be  taken  into  account. 

When  developing  a  solution  method  in  more  than  a  single  direction  it  is  often  useful 
to  define  a  localized  coordinate  system  for  each  discrete  point  in  the  domain  |2l].  The 
conservative  variables  and  fluxes  can  be  calculated  as  a  function  of  the  unit  vectors  of  the 
local  coordinates.  For  an  arbitrary  curvilinear  domain  the  local  normal  coordinate  system  at 
each  point  or  control  volume  will  not  necessarily  align  itself  with  the  cartesian  coordinates 
system.  However,  a  coordinate  transformation  can  be  developed  such  that  the  physical 
domain  is  mapped  to  the  local  computational  (cartesian)  domain.  This  transformation  is 
done  through  the  use  of  “grid  metrics” ,  defined  as 

tx  =  dy 

J  drj 
_  _  dx 

J  drj 

rh=_dy  (2.9) 

J 

rjy  _dx 

7 

J  =CxVy  £,yVx 
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where  the  \x,  y\  domain  is  the  physical  domain  and  the  [£,  p]  is  the  computational  domain. 
The  term  J  in  two  dimensions  can  be  considered  as  the  inverse  of  the  cell  area.  It  is  often 
convenient  to  specify  that  the  computational  domain  is  uniform  such  that  =  A?/  = 
1.  In  this  manner  a  non-uniform  curvilinear  physical  domain  can  be  transformed  into 
a  uniform  rectilinear  computational  domain.  Using  these  definitions  for  grid  metrics  the 
Euler  equations  can  be  written  in  generalized  coordinate  differential  form  as 


1  dU  8F  dQ 
J  dt  dp 


(2.10) 


where  T  and  Q  are  the  conservative  flux  vectors  such  that 


F  =  ^jF  +  ^!-G 

U  J 

g=rhF+rhG  (2.H) 

J  J 


The  conservative  variable  vector  U  and  flux  vectors  F  and  G  can  also  be  expressed  in  terms 
of  the  primitive  variables 


p 

pu 

pu 

pu 

,F  = 

pu2  +  p 

,G  = 

puv 

pv 

puv 

pv2  +  p 

_  Ft  _ 

_  ( Et+p)u  _ 

_  (Et  +  p)v  _ 

(2.12) 


Finally,  the  equation  of  state  is  again  used  for  closure  such  that  with  an  ideal  gas  total 
energy  may  be  calculated  as 


Et  = 


P 

7-1 


+  2  p(«2  +  «2) 


(2.13) 


As  with  one-dinrension,  it  is  often  beneficial  to  express  the  Euler  set  in  conservative  integral 
form  for  use  in  a  finite-volume  method.  Using  the  same  definition  of  p.  u,  v  and  p  it  can  be 
shown  that  the  corresponding  integral  form  is 


d_ 

dt 


UcN  + 


J (Fi  +  Gj)  ■  dS  =  0 


(2.14) 
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Pi  ; 

p2 

p'  = 

P2  [— > 

VTi  ’ 

Vni  = 

VN2*0 

(a)  Contact  Surface  (b)  Slip  Line  (c)  Shock  Surface 

Figure  2.1:  Discontinuous  Solutions  to  the  Rankine-Hugoniot  Relations 

Eulers  equations  allow  for  both  strong  (continuous)  and  weak  (discontinuous)  solutions 
|11|.  Strictly  speaking,  the  differential  form  of  Eulers  equations  requires  a  solution  with 
smoothness  or  one  with  continuous  derivatives.  In  order  to  capture  the  weak  solutions 
an  integral  form  must  be  used.  It  can  be  shown  that  the  Euler  set  allows  three  types  of 
discontinuous  solutions,  all  of  which  satisfy  the  Rankine-Hugoniot  jump  conditions:  contact 
surfaces,  slip  lines,  and  shock  surfaces  (Figure[2.1[).  Contact  surfaces  are  defined  where  there 
is  a  discontinuity  in  density,  but  pressure  and  velocity  are  continuous.  At  the  interface 
between  the  two  states  a  boundary  exists  with  zero  mass  flow  across  it.  Slip  lines  are 
similar  to  contact  surfaces,  but  there  also  exists  a  discontinuity  in  tangential  velocity  across 
the  boundary  while  pressure  and  normal  velocity  remain  continuous.  The  third  solution, 
with  mass  flow  across  the  boundary,  is  a  shock  surface  where  pressure,  density,  and  normal 
velocity  all  experience  discontinuities.  For  the  shock  surface  tangential  velocity  remains 
continuous  across  the  shock. 

Since  stagnation  pressure  is  not  constant  across  a  shock,  entropy  likewise  is  discontin¬ 
uous.  Mathematically  the  Euler  equations  allow  for  the  existence  of  a  condition  where  the 
change  in  entropy  across  a  discontinuity  could  be  negative,  forming  an  “expansion  shock” 
of  sorts.  However,  this  result  would  be  in  violation  of  the  second  law  of  thermodynamics: 
that  the  change  in  entropy  of  a  system  must  always  be  greater  than  or  equal  to  zero.  A 
valid  solution  technique  must  therefore  be  able  to  enforce  this  constraint  due  to  the  second 
law  thus  avoiding  non-physical  behavior. 

In  the  next  chapter  the  numerical  schemes  to  solve  the  linear  wave  equation  and  the 
Euler  equations  will  be  presented. 
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III.  Numerical  Methods 


The  Linear  Wave  Equation  and  Euler’s  Equations  (both  one  and  two  dimensional)  are 
hyperbolic  PDE’s  of  two  or  more  independent  variables.  In  order  to  reach  a  numerical 
solution  they  must  be  separated  into  two  distinct  portions:  spatial  terms  (x,  y )  and  temporal 
terms  (t).  A  common  solution  method  consists  of  first  calculating  a  linearized  spatial  term 
at  a  given  time  then  integrating  explicitly  with  respect  to  time  to  reach  a  solution  at  the 
next  time  step.  This  process  is  then  repeated  until  the  desired  end  state  is  reached.  In 
order  to  calculate  the  spatial  terms  the  equations  must  first  be  discretized  or  simplified  for 
evaluation  at  a  finite  number  of  discrete  points  in  a  computational  domain. 


3.1  Formulation  of  the  Compact  Difference 

The  compact  difference  begins  with  a  generalized  7-point  implicit  stencil  1 12 1 : 


nr/  |  /■/  .  rl  .  rl  ,  n  rl  ./)  ‘  1  fi — 1  .  ifi-\-2  fi  —  2  .  f i-\-  3  fi  —  3  /q  i\ 

Pfi-2  +  a/i-i  +  fi  +  afi+ 1  +  Pfi+2  =  a - ^ - +  b - — - +  c - — -  (3.1) 


AAx 


6Ax 


where 


,/  _  dfj 

Ji  dx 


(3.2) 


If  TSE’s  for  f,  fi+ 1,  etc.  are  substituted  into  Equation  3.1  a  system  of  equations  can 
then  be  written  such  that  the  coefficients  a,b  and  c  cause  desired  terms  of  the  TSE  to 
cancel.  Solving  this  system  of  equations  will  produce  the  coefficients  needed  for  an  implicit 
expression  for  f  with  a  truncation  error  consisting  of  the  remaining  terms.  Using  this 
method  the  fourth  order  representation  of  the  first  derivative  becomes: 


I  rl  |  rl  |  I  rl  ^  1  fi — 1  .  A  \4 

\fi~  i  +  ^  +  I^+i  =  4 - ~Ax - +  °(AX) 


(3.3) 


This  particular  compact  stencil  will  be  henceforth  referred  to  as  the  “cC”  compact  scheme. 
Note  that  the  stencil  has  been  reduced  from  5  points  for  a  traditional  finite  difference  (i- 


2  through  i+2,  see  Equation  1.3)  to  3  points  (i-1  through  i+1).  Other  values  for  these 


coefficients  and  their  respective  orders  are  summarized  in  Table  |3.1[ 

These  formulations  are  all  based  on  the  assumption  that  the  value  of  the  function  / 
and  the  desired  derivative  f  reside  at  the  same  discretized  locations.  This  methodology  can 


3-1 


also  be  used  when  /  and  f  are  not  co-located  assuming  the  appropriate  TSE’s  are  used.  For 
example,  in  a  finite-volume  discretization  scheme  the  values  of  a  flux  /  might  be  known  at 
cell  interfaces  while  the  change  in  flux  df/dx  would  be  desired  at  the  cell  centers.  TSE’s  for 
the  right  hand  side  would  be  expanded  about  i  ±  1/2  whereas  TSE’s  for  the  left-hand  side 
would  be  expanded  about  i  —  1,  i,  and  i  +  1.  Again,  a  system  of  equations  can  be  written  in 
terms  of  a,  (3,  a,  b,  and  c  such  that  desired  terms  in  each  TSE  cancel.  The  remaining  TSE 
terms  make  up  the  truncation  error.  Using  this  method  a  fourth  order  set  of  coefficients 
using  these  interface  values  could  be  calculated  using: 


af'_1  +  f'  +  Pf'+1  =  —f, 


'i+i 


Ax^1  2 


(3.4) 


with 


3 

a  =-(a  +  P) 

b  =\{a  +  (3) 

o 


All  of  these  stencils  produce  an  implicit  formulation  for  the  value  of  f .  The  resulting 
solution  matrix  is  sparse,  but  banded.  The  left-hand  side  coefficients  a  and  [3  are  free 
parameters  which  determine  the  solution  matrix  type.  For  (3  =  0  and  2a  <  1  the  resulting 
matrix  is  diagonally  dominant  and  tri-diagonal  which  can  be  easily  solved  using  the  Thomas 
algorithm  1 10  .  For  eighth  and  higher  order  schemes  (3  ^  0  and  the  resulting  matrix  is  penta- 
diagonal  which  must  be  solved  using  other  methods. 

For  this  thesis,  the  fourth-order  version  of  Equation [3T] was  used  for  the  cC  scheme  and 
the  implicit  operator  was  a  diagonally-dominant  tri-diagonal  system.  A  low-storage  version 
of  the  Thomas  algorithm  pi]  was  then  employed.  The  Thomas  algorithm  is  essentially 


Table  3.1: 


Various  Values  of  Parameters  from  Equation  (3.1) 


a 

P 

a 

b 

c 

Order 

a 

0 

§(a  +  2) 

5  (4a  —  1) 

0 

4th 

a 

0 

g(a  +  9) 

15  (32a -9) 

ro(!-3a) 

6th 

a 

'oo 

e 

I 

03 

|(12 -7a) 

jgg(568a  -  183) 

5o(9a  —  4) 

8th 
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Gaussian  elimination  modified  to  take  advantage  of  the  sparse  matrix.  The  implicit  operator 
is  first  put  into  upper  triangular  form  by  computing  new  diagonal  elements.  The  second 
step  computes  new  right-hand  values  based  on  this  upper  triangular  matrix.  The  final  step 
is  to  back-substitute  and  compute  the  unknown  values.  Since  the  only  non-zero  terms  of  the 
left-hand  matrix  are  the  lower  diagonal,  diagonal,  and  upper  diagonals,  they  are  the  only 
terms  passed  to  the  Thomas  subroutine  in  the  form  of  three  column  vectors.  The  Thomas 
algorithm  has  the  advantage  of  solving  tri-diagonal  system  very  quickly  with  a  minimal 
number  of  calculations  and  storage. 


3.2  Boundary  Conditions 


Compact  methods  such  as  Equations  |3.1|  and  3.3  use  a  central-difference-like  stencil 


which  can  lead  to  difficulties  near  boundaries.  For  example,  Equation  3.3  is  only  usable  in 
the  domain  from  2  <  i  <  n  —  1  where  n  is  the  number  of  known  values  of  /.  At  the  points 
i  =  1  and  i  =  n  another  method  must  be  used  since  the  stencil  would  reach  outside  of  the 
solution  domain. 


For  Dirschlet  and  Neumann  boundary  conditions  ghost  cells  can  be  readily  extrapo¬ 
lated  from  the  domain  and  used  as  known  points.  For  these  cases  the  values  of  /  or  f  would 
be  readily  known  and  the  value  could  be  moved  to  the  right  hand  side  of  the  equation. 


An  alternate  method  is  to  use  one-sided  differences  which  are  derived  in  a  similar 


manner  as  Equation  3.1  In  this  case,  the  left  hand  stencil  does  not  extend  out  of  the 
known  domain  (1  <  i  <  ni).  A  one-sided  approximation  to  f  at  a  minimum  boundary 
(i  =  1)  could  be  used  as: 


f'i  +  a/i+l  =  +  bfi+l  +  °fi+ 2  +  dfi+ 3  (3-5) 

With  [a,  a,b,  c,  d]  =  [-17/6,  3/2,  3/2,  -1/6]  the  solution  is  fourth-order.  It  is  worth  noting 
that  when  the  solution  becomes  one-sided  the  stencil  width  of  the  right  hand  side  increased 
from  three  points  to  four.  A  tradeoff  to  widening  the  right  hand  side  would  be  to  reduce 
the  order  of  accuracy  of  the  approximation  at  the  boundary. 
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3.3  Application  to  the  ID  Linear  Wave  Equation 


As  previously  stated  the  Linear  Wave  Equation  is  given  by: 


df  ,  df_ 
~di  +  ad^~° 


(3.6) 


One  possible  explicit  solution  to  the  linear  wave  equation  is  the  first-upwind  formulation  1 12 1 : 


fn+ 1  fn  a  At  i  j,n  ^ 

n  ~  Ji  Ji~l) 


(3.7) 


Known  values  of  /  at  time  n  are  used  to  explicitly  calculate  new  values  of  /  at  n+1.  In  effect 
both  temporal  and  spatial  time  steps  are  calculated  at  the  same  time.  This  method  is  first 
order  accurate  in  time  and  space  (0(Ax,  At)).  While  this  method  is  simple  to  implement 
it  introduces  a  great  deal  of  artificial  dissipation  for  values  of  v  <  1  (where  v  =  aAt/Ax) 
which  results  in  a  diminishing  wave.  A  more  accurate  solution  would  result  from  using  a 


compact  difference  method  (Equation.  3.3)  for  the  spatial  term  df  /dx  and  then  a  separate 
time  integration  method  for  the  temporal  term.  To  make  the  spatial  integration  fourth-order 
the  following  formulation  was  used: 


rt  .  rf  i  w  fi+1  fi — 1  .  t  ./f- 1~2  fi — 2 

afi-1  +  fi  +  afi+ 1  =  a - — - h  b- 


2Ax 


4Ax 


(3.8) 


where 


a  =3^a  +  2^ 

b  =-(3a  —  1) 
4  7 

1 

a  =  — 

10 


(3.9) 


Note  that  Equation  3.9  is  used  instead  of  Equation  3.1  so  that  the  left-hand  side  forms  a 
tridiagonal  system  of  equations.  Known  values  of  /  at  time  n  are  used  to  make  the  right 


hand  side  of  Equation  3.9  then  the  implicit  problem  was  solved  for  f  using  the  Thomas 
algorithm  1 21  .  The  new  value  of  /  at  time  n  +  1  was  then  determined  using  a  corresponding 
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fourth-order  Runge-Kutta  time  integration  method: 


h=r 
h  =r  - 
h  =r  - 
u  =r  - 


a  At  df\ 

4  dx 
aAt  3/2 
3  dx 
aAt  dh 
2  dx 


(3.10) 


fn+1  _  aAtdh 

dx 

making  the  solution  globally  fourth-order  (0(Ax4,  At4)).  Note  that  the  implicit  system 


formed  by  Equation  3.9  needs  to  be  solved  at  each  stage  of  the  Runge-Kutta  scheme. 


3-4  Application  to  the  ID  Euler  Equations 


The  flux  term  of  the  one-dinrensional  Euler  set  (Equation  2.4)  is  a  nonlinear  trans- 
portive  term  similar  to  the  df  /dx  term  of  the  LWE.  The  temporal  term  can  be  calculated 
using  the  Runge-Kutta  time  integration  scheme  as  was  used  on  the  LWE,  but  the  flux  term 
cannot  be  solved  using  the  same  compact  difference  method  due  to  the  inherent  design  of 
the  compact  difference  stencil.  A  major  disadvantage  of  the  formulations  such  as  Equation 


3.1  or  equation  3.3  is  that  they  include  central-like  difference  terms  on  the  right-hand  side  of 
the  formulation.  Central  differences  when  used  on  nonlinear  hyperbolic  problems  produce 
unconditionally  unstable  oscillations  around  discontinuities.  This  was  not  an  issue  with  the 
LWE  which  was  a  scalar  linear  problem  with  a  single  characteristic  of  a  constant  sign.  The 
Euler  set  in  one  dimension  has  three  characteristics:  u,  u  +  a,  and  u  —  a,  (where  a  is  the 
speed  of  sound).  The  characteristics  are  the  eigenvalues  of  the  flux  Jacobian  and  dictate 
the  direction  and  speed  of  information  propagation.  For  a  subsonic  problem  two  of  the 
characteristics  will  propagate  in  the  flow  direction  while  one  will  propagate  in  the  opposite 
direction  of  the  flow.  For  a  supersonic  problem  all  three  characteristics  run  in  the  flow 
direction.  To  maintain  stability  the  solution  must  take  the  direction  of  these  characteristics 
into  account  and  only  use  values  of  the  function  in  the  upwind  or  opposite  direction  of  the 
characteristic. 
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3.5  Roe’s  Approximate  Riemann  Solver 

To  develop  an  upwind  solution  method  the  one-dimensional  Euler  equations  can  be 
rewritten  in  non-conservative  quasi-linear  form  in  terms  of  the  Jacobian  A: 


dU 

~di 


AdU  n 
+  —  0 
ox 


(3.11) 


where  the  flux  Jacobian  is 


.  dF 
~  dU 


(3.12) 


Philip  Roe  11  postulated  that  if  A  were  evaluated  at  some  average  interface  value  (A) 


rather  than  at  cell  centers  a  conservative  upwind  scheme  could  be  developed  such  that: 


Fi+i=-(Fi  +  Fi+1-\A\(Ui+1  -170) 


(3.13) 


Roe  noted  that  both  the  conservative  variable  vector  U  and  the  Euler  flux  vector  F  were 
quadratic  functions  of  a  parameter  vector  z  such  that 


z  =  Vp 


u 


H 


(3.14) 


where  p  is  density,  u  is  x-velocity,  and  H  is  total  enthalpy.  Roe  then  showed  using  the 
properties  of  quadratic  functions  that  the  Roe  matrix  A  is  exactly  the  flux  Jacobian  A, 
evaluated  with  Roe-averaged  variables.  For  a  one-dinrensional  system  Roe’s  scheme  begins 
with  “Roe-averaged”  primitive  variables  being  calculated  at  cell  interfaces  using  the  square 
root  of  the  product  of  the  densities  R 


Ri+l/2 


Pi+ 1/2 
&1+1/2 

Hi+l/2 


=Ri+l/2Pi 

Ri+l/2ui+l  +  ui 

Ri+ 1/2  +  1 

-^i+l/2-^i+l  +  Hi 
Hi+l/2  +  1 


(3.15) 
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The  average  speed  of  sound  at  the  interface  is  then  calculated  from  these  values  with 


a  =( 7-1)  (tf-y 

The  associated  eigenvalues  and  right  eigenvectors  to  the  Roe  matrix  A  are 

Am  =  u  A^2)  =  u  T  a  Am  —  ^  —  u 

and 


(3.16) 


(3.17) 


f^)  = 


1 

1 

1 

u 

~  2 

,  f(2)  = 

u  +  a 

II 

CO 

u  —  a 

U 

2 

FI  +  ua 

H  —  ua 

P_ 
2  a 


(3.18) 


Next,  the  wave  amplitudes  at  the  interface  are  calculated  using  simple  averages 


5w\  = 

dp 

dp  -  ^ 

az 

dw  2  = 

dp 

ou - 

pa 

5u’3  = 

dp 

ou - 

pa 

dui+ 1/2  = 

^i+l 

dPi+ 1/2  = 

Pi+l  ~  Pi 

dPi+ 1/2  = 

Pi+ 1  Pi 

(3.19) 


and  finally,  the  numerical  flux  F  is  found  from 


Fi+ 1/2 


(Fi  +  Fi+ 1)  -  ^ 


(3.20) 


where  Fi  and  Nj+i  are  the  exact  fluxes  based  on  the  left  and  right  cell  states.  When 


written  in  the  form  of  Equation  3.20  the  summation  term  can  be  thought  of  as  an  artificial 
dissipation-like  term  which  was  calculated  using  the  appropriate  upwind  information. 

One  drawback  to  Roe’s  scheme  as  presented  above  is  that  the  solution  does  not  see 
the  sonic  point  and  will  admit  an  expansion  shock  as  a  valid  solution  [11  .  This  non-physical 
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behavior  can  be  eliminated  by  the  addition  of  slight  artifical  dissipation  in  the  calculation 
of  the  A  eigenvalues.  Harten  and  Hyman  |8  suggest  the  following  entropy  correction: 


I  A  |  mod  * 


1  'Mi+1/2 

1 

2 


\2 

i  + 1/2 


^  l^li+1/2  >  e 
if  l^li+1/2  <  e 


where  e  is  a  value  such  that 


(3.21) 


e  =  max  [0,  (Ai+1/2,  -A*),  (Am  -  Ai+i/2)] 


(3.22) 


3.6  Roe’s  Scheme  in  Two  Dimensions 

When  Roe’s  solution  is  extended  to  two  dimensions  the  procedure  needs  to  be  modified 
slightly  (23  .  Since  there  is  more  than  one  dimension  the  contravariant  velocities  U  and  V 
must  be  used  instead  of  the  vector  quantities  u  and  v. 


U  =Z,xu  +  iyV 

V  =rjxu  +  rjyV 


m\ 

t  = 

'y  |V£| 


Vx 


Vx 

|Vr7| 


Vy 


_  Vy 
I  V?7| 


(3.23) 


where  U  and  V  are  the  velocity  components  in  the  £  and  r]  directions,  respectively.  The 
flux  vector  in  the  £  direction  is 


7  = 


M 

j 


pU 

pUu  +  £xp 

pU'V  +  iyP 

(■ Et+p)U 


(3.24) 
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The  Roe  vector  which  is  a  product  of  the  Roe  matrix  and  the  change  of  conserved  variables 
across  the  cell  interface  is  given  by 


\A\(ui+l  -  Ui))  = 


|V£|  ua4  +  j^a5  +  a6 
J  va4  +  |^ya5  +  a7 


(3.25) 


Ha4  +  Ua 5  +  ua  6  +  va7  —  (-y-r 


where 


=\U\{^P-% 

az 

012  =  2^  +  “I  i^P  +  paAU) 

«3  =  2^2  1^  “  “I  (Ap  ~  paAU) 

O4  —0(1  +  02  +  0(3 
a  5  =a(a2  —  a  3) 


(3.26) 


«6  =\U\{pAu  -  ^-pAU) 
a7  =\U\(pAv  -  T^i  pAU) 

U  - £x  U  ^yV 


Similarly,  the  flux  vector  in  the  p  direction  is 


|  Vr?|  pVu  +  rjxp 

J  pVv  +  pyp 

( Et+p)V 


(3.27) 


The  Roe  matrix  in  this  direction  is  given  by 


ia‘  +  i«r“5  +  “« 

J  va4  +  iw"5  +  °7 


(3.28) 


Ha  4  +  Vos  +  -ua6  +  iia7  —  (wb[)ai 
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«!=|V|(A  p-^4) 

az 

°2  =2#  ^  +  a\(Ap  +  paAV) 
«3  =^2 1^  ~  o|(Ap  -  paAV) 

Ct 4  =OL\  +  012  H“  013 
«5  =a(oi2  —  OL  3) 

a6  =|V|(/5Au  -  j^|pAV) 
a7  =|V|(/5Av  -  pAV ) 

V  =pxu  +  PyV 


(3.29) 


The  Roe  averaged  variables  are  defined  in  a  similar  manner  as  the  one-dimension  case, 


R 


i+1/2 


Pi+ 1 
Pi 


Pi+ 1/2  =Ri+l/2Pi 

Ri+l/2ui+l  +  ui 

Ri+ 1/2  +  1 

-^i+l/2vi+l  +  vi 
^i+1/2  +  1 
-^i+l/2-^i+l  +  Hi 

R-i+1/2  +  1 


Mi+l/2 

Vi+l/2 


(3.30) 


H, 


i+1/2 


Finally,  the  Roe  fluxes  are  calculated  using 


^+i  =  J(^  +  ^+i  -  |i|(t/m  -  I7f))  (3-31) 

2  A 

Qj+l=l-(Qj  +  Qj+i-\A\(Uj+1-Uj)) 


The  same  entropy  fix  as  stated  in  Equation  3.21  was  also  used  to  eliminate  the  possibility 
of  non-physical  expansion  shocks  in  the  2D  Roe  flux  calculation. 
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3. 7  Higher-  Order  Approximations 

Roe’s  method  as  presented  above  is  a  first  order  approximation  in  space.  In  practice 
this  method  is  too  dissipative  to  be  used.  Odd-ordered  methods  tend  to  have  dissipative 
qualities  to  them  due  to  the  dominance  of  even-ordered  derivative  terms  in  the  truncation 
errors.  This,  combined  with  the  fact  that  Roes  method  uses  an  artificial  dissipation-like 
term  to  calculate  the  upwind  fluxes,  results  in  a  smearing  of  gradients  and  a  loss  of  signal 
magnitude  in  solutions. 

An  improvement  comes  through  the  use  of  Roe’s  method  in  conjunction  with  a  variable 
extrapolation.  Roe’s  method  is  essentially  a  weighted  average  of  the  flux  at  a  cell  interface 
expressed  as  a  function  of  the  primitive  variable  vector  V  ([p,  u,p]T)  in  each  adjoining  cell. 
A  first  order  approximation  assumes  that  the  value  of  the  primitive  variable  is  constant 
throughout  the  entire  cell  with  variation  only  at  the  cell  interfaces.  If  variation  within  the 
cell  is  allowed  a  higher  order  reconstruction  can  be  developed  to  obtain  more  accurate  left 
and  right  primitive  variable  states  at  the  cell  interface. 

VanLeer’s  Montone  Upwind-centered  Scheme  for  Conservation  Laws  j;22|  (MUSCL) 
defines  a  series  of  piecewise  linear  approximations  to  the  solution  between  two  cells  at  the 
interface.  Instituting  a  TSE  about  the  cell  average  values  for  a  primitive  variable  v  yields: 

vi  =  vi  +  (x-xi)  —  \i  +  ^2 -  -^It  +  OiAx3)  (3.32) 


where 


Vi  =  -£-  [Xi+h  vdx 
Ax  J: r.  ! 


(3.33) 


Upon  simplification,  the  value  for  v  can  be  written  as 


_  .  dv  1  d2v 

Vi  =  Xi  +  (x-Xi)-  +  -  — 


(x  -  Xi )2  - 


Ax2 


(3.34) 


Van  Leer  chose  to  scale  the  last  term  of  this  approximation  to  allow  a  family  of  piecewise 
linear  and/or  quadratic  extrapolation.  If  central  difference  approximations  are  used  for  the 
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first  and  second  derivatives  the  left  state  cell  interface  value  of  v  can  be  calculated  using 


vl  =  Vi  +  ^(1  -  n)(vi  -  Vi- 1)  +  |(1  +  n)(vi+ 1  -  Vi) 


(3.35) 


and  the  right  state  cell  interface  value  of  v  using 

vR  =  vi+i  -  |(1  +  n)(vi+ 1  -  Vi)  -  ^(1  -  K)(vi+2  -  fj+i)  (3.36) 


The  choice  of  k  determines  the  order  of  the  extrapolation,  k  =  1/3  leads  to  a  third-order 
accurate  extrapolation,  whereas  n  =  —  1  would  lead  to  a  3  point  fully  upwind  discretization. 


To  make  Equation  3.32  higher  order,  substitute  L  and  R  states  for  i  and  i  +  1  respectively. 


For  this  thesis,  the  third-order  k  =  1/3  method  will  be  used  and  henceforth  referred  to  as 
“oRoe” ,  or  “Original  Roe” . 


3.8  Variable  Limiters 

Total  variation,  or  the  total  change  in  any  general  solution  /  is  defined  as 


TV  = 


df 


dx 


dx 


(3.37) 


If  this  total  variation  is  bounded  in  both  time  and  space,  such  that 


TV(/n+1)  ^  TV(^ 


(3.38) 


the  solution  is  said  to  be  total  variation  diminishing,  (TVD)  and  the  following  properties  of 
nronotonicity  hold  true  [II]:  TVD 

1.  No  new  local  extrema  can  be  created 

2.  The  value  of  a  local  minimum  is  non-decreasing  and  the  value  of  a  local  maximum  is 
non-increasing 

If  these  conditions  are  met  the  scheme  is  considered  monotone.  In  essence  a  TVD  scheme  is 
stable  such  that  an  introduced  perturbation  (due  to  error,  etc)  in  the  solution  will  not  grow 
without  bounds,  and  the  solution  will  remain  bounded.  While  the  use  of  variable  extrap- 
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olation  is  useful  for  increasing  the  order  accuracy  of  a  solution  it  can  introduce  additional 
errors  due  to  aliasing  and  Gibb’s  effects  at  gradients  and  discontinuities  [l8j.  These  errors 
cause  solutions  to  become  non-TVD  and  often  cause  the  solution  to  diverge.  For  this  reason 
a  limiter  is  often  introduced  to  enforce  the  TVD  property. 

There  are  various  limiter  functions  which  all  have  a  similar  purpose:  reduce  or  elimi¬ 
nate  non-physical  oscillations  near  solution  extrema.  In  order  to  work  properly  the  limiter 
must  be  used  in  conjunction  with  a  scheme  that  is  TVD  when  first  order.  The  limiter  func¬ 
tions  by  extending  the  scheme  to  a  higher  order  when  the  solution  is  smooth  then  reverting 
back  to  the  first-order  TVD  when  the  solution  is  non-smooth.  One  limiter  function  which 
is  commonly  used  is  the  minmod  limiter,  defined  as  (Tl] 


minmod(x',  y) 


x  if  |x|  <  \y\  and  xy  >  0 
<  y  if  |x|  >  \y\  and  xy  >  0 
0  if  xy  <  0 


In  compact  notation  this  limiter  is  often  expressed  as 


(3.39) 


minmod(x,y)  =  sign(x)  •  nrax[0,  min(|x|, sign(x  •  y)\ 


(3.40) 


In  effect  the  limiter  is  a  switching  function  to  determine  how  the  variable  extrapolation  is 
to  be  constructed.  Inherent  to  their  design  limiters  have  the  effect  of  reducing  the  order 
of  accuracy  of  the  solution  near  extrema.  Because  of  this  behavior  variable  extrapolation 
schemes  are  not  of  a  global  order  accuracy.  This  degradation  of  order  is  often  manifest  as 
a  broadening  for  a  shock  or  a  smearing  for  a  contact  surface. 


3.9  Upwind  Compact  Formulations 

If  an  upwind  bias  could  be  incorporated  into  a  compact-difference  method  a  solution 
with  global  higher-order  accuracy  could  theoretically  be  developed.  If  the  upwinding  were 
done  in  the  proper  manner,  non-physical  behavior  could  be  avoided  and  the  use  of  limiters  or 
filters  would  be  unnecessary.  As  previously  mentioned  it  was  suggeted  by  David  Gottlieb  [7] 


that  a  combination  of  a  compact  scheme  like  Equation  3.1  and  Roe’s  scheme  could  create  a 
higher-order  solution  with  minimal  effort. 
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An  investigation  of  this  theory  shows  that  it  is  only  partially  correct.  Deng  and 
Maekawa  |3  show  this  by  beginning  with  the  definition  of  Roe’s  scheme, 

£+1/2  =  \[F{Ur)  +  F(UL)  -  \A\(UR  -  UL)]i+ 1/2  (3.41) 

Assuming  the  left  and  right  states  Ul  and  Ur  are  interpolated  to  rth-order  accuracy  the 
Roe  interface  flux  becomes 


h+i/2  =  F(Ui+ 1/2)  +  d(xi+1/2)hr  +  0(hr+1)  (3.42) 


and  the  characteristics  of  the  variable  extrapolation  scheme  and  the  cell-centered  differenti¬ 


ation  scheme  become  one.  If  Equation  3.42  is  substituted  back  into  the  compact  formulation 


of  Equation  1.5  and  TSE’s  applied  at  each  point  i,  an  expression  for  F'  can  be  shown  to  be 


.  d(xi+1/2)-d(xi_1/2)  17  (d5P\ 

+  “ - h - k  +  5280  (  W>) 


h4  +  0(h6,hr+1 ) 


(3.43) 


Assuming  that  d(x)  is  continuous  the  value  of  F'  will  be  fourth  order  accurate  if  r  >  4. 
If  r  <  4  the  value  of  F'  will  be  reorder  accurate.  Hence,  if  the  Roe  fluxes  are  only 
calculated  as  first  order  the  resulting  compact  scheme  will  produce  a  first  order  solutions.  If 
variable  extrapolation  (MUSCL)  is  used  to  increase  the  accuracy  of  the  left  and  right  states 
to  third-order  the  resulting  F’  will  also  be  thrid-order.  It  is  worth  noting,  however,  that 
Deng  and  Maekawa  |3  did  show  that  the  once  unconditionally  unstable  method  of  Equation 
P  was  now  stable  without  the  use  of  limiters  or  filters.  This  is  due  to  the  fact  that  the 
Roe  flux  calculations  are  taking  the  characteristics  into  account  and  are  in  effect  making  a 
compact-upwind  method. 

The  next  logical  step  is  to  try  and  increase  the  order  of  the  variable  extrapolations  at 
the  left  and  right  interface  states  to  at  least  fourth-order.  This  can  be  accomplished  through 
the  use  of  an  implicit  interpolation  similar  in  form  to  a  compact  difference.  However,  if  the 
stencil  of  the  interpolant  includes  a  discontinuity  the  solution  will  include  nonphysical  oscil¬ 
lations.  If  the  stencil  could  be  dynamically  chosen  such  that  the  interpolant  stencil  avoids 
discontinuities  higher-order  reconstruction  can  be  developed  without  the  use  of  filtering  or 
limiting. 
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3.10  Implicit  Variable  Interpolation 

Deng  and  Maekawa  |3|  developed  a  method  by  determining  the  smoothness  of  the 
function  that  is  being  interpolated.  The  absolutes  of  the  first  and  second  Essentially  Non- 
oscillatory  (ENO)  differences, 


Dij  —  vj+ 1  Vj  j 

D2j  =  Vj+ 1  —  2  Vj  +  Vj- 1 


(3.44) 


can  be  used  as  a  measure  of  the  smoothness  of  the  stencil  |9|,  where  v  is  any  primitive 
variable.  The  first  ENO  difference  is  akin  to  the  slope  of  the  function  while  the  second  ENO 
difference  is  the  curvature  of  the  function.  The  smaller  \Di\  and  \D2\  the  smoother  the 
function.  Using  these  definitions  along  with  a  weighting  factor  k,  a  switch  function  K  for 
the  left  and  right  states  can  be  developed  so  the  compact  stencil  can  be  adjusted  to  avoid 
discontinuities  [3  .  For  the  left  state: 


Kl 


1,  if  kIDij-iI  <  l-Dijj,  \D2j_i\  <  n\D2j\\ 

<  3,  if  \D\j  <  \D2jjri\  <  n\D2j\\ 

2,  otherwise. 


(3.45) 


The  interpolation  for  vLj+1/2  is  defined  as 


avLj-i/2  +  VLj+i/2  +  PvLj+3/2  =  avj-3  +  bvj-2  +  cvj-1  +  dvj  +  evj+i  +  fvj+2  +  gvj+3  (3.46) 

The  resulting  values  for  a,b,c,d,e  for  a  fourth  order  interpolate  of  vl  are  summarized  in 
Table  lT2l 

Similarly  for  the  right  state,  the  interpolation  for  vrj+ \/2  is  defined  as 


avRj-1/2  +  vRj+i/2  +  f3vRj+3/2  =  avj-3  +  bvj-2+cvj-l  +  dvj+evj+i  +  fvj+2+gvj+3  (3.47) 
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Table  3.2:  Various  Values  of  Parameters  from  Equation  (3.46) 


Kl 

a 

p 

a 

b 

c 

d 

e 

f 

g 

1 

a 

p 

T(a-  5) 

T(21  -  5a) 

T  (15a -35) 

(5  a  +  35) 

0 

0 

0 

2 

1 

2 

l 

10 

0 

0 

1 

10 

1 

1 

2 

0 

0 

3 

a 

p 

0 

0 

0 

Tg(S-P) 

Yg(!5  +  9/3) 

33(90-5) 

M'-P) 

Table  3.3: 


Various  Values  of  Parameters  from  Equation  (3.47) 


Kr 

a 

p 

a 

b 

c 

d 

e 

f 

g 

1 

a 

p 

0 

0 

0 

T  (5/3  +  35) 

33(15/3  —  35) 

33(2! -5 P) 

Ie(P  -5) 

2 

l 

10 

1 

2 

0 

0 

1 

2 

1 

1 

10 

0 

0 

3 

a 

P 

TgO-a) 

Tg(9a-  5) 

T(15  +  9a) 

33(5  -  a) 

0 

0 

0 

with  the  switch  function  Kr  defined  as 


Kr 


1,  if  n\Dij+i\  <  \Dij\,  | |  <  R\D2j+i\', 

<  3,  if  \D\j  <  n\D\j+i,  \D2j\  <  n\D2j+i\-, 

2,  otherwise. 


(3.48) 


and  the  corresponding  coefficients  determined  by  Kr  shown  in  Table  3.10 


The  parameters  a  and  f3  are  free  parameters  which  determine  the  form  of  the  left- 
hand  side  of  the  implicit  formulation.  For  this  work  they  were  chosen  such  that  a  +  f3  <  1 
to  maintain  a  diagonally  dominant  system  that  could  be  also  solved  with  the  Thomas 
algorithm  |2l] .  Deng  and  Maekawa  |3  suggest  the  weighting  factor  k  should  be  chosen 
0  <  k  <  1,  such  that  Kr  =  2  and  Kr  =  2  (a  central-type  interpolation)  is  used  as  much  as 
possible.  Since  in  essence  this  weighting  factor  determines  when  the  solution  is  sufficnetly 
non-snrooth  to  warrant  a  one-sided  interpolant  it  is  anticipated  that  the  optimal  value  of  k 
will  vary  with  each  problem  the  scheme  is  applied  to. 


This  scheme  is  called  a  “Compact  Adaptive  Interpolation  Scheme”  by  Deng  and 
Maekawa  |3  .  For  the  purpose  of  this  thesis,  it  will  be  referred  to  as  “cRoe”,  or  the  “Compact 
Roe”  scheme. 
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3.11  Time  Integration 


Once  the  numerical  Roe  flux  is  known  a  finite-volume  scheme  was  employed  such  that 
the  residual  R  is  defined  as 

Ri  =  Fi+ 1/2  —  Ft- 1/2  (3.49) 

in  one  dimension,  or 


Ri,j  — 


F i+l/2,j  -  F j— l/2,j  ftj+1/2  -  Qi,j- 1/2 


A? 


At/ 


(3.50) 


in  two  dimensions.  For  the  compact  schemes,  solving  the  implicit  matrix  results  in  deriva¬ 
tives  of  the  fluxes  directly  so  that  the  residual  becomes 


R  =d^+9l 

1  d£  dr] 


(3.51) 


Once  the  residual  is  determined  the  time  integration  of  the  dU/dt  term  can  then  be  ac¬ 
complished  using  a  four-stage  CFD-type  Runge-Kutta  scheme  jl3  .  This  method  explicitly 
solves  for  the  conserved  variables  at  time  level  n  +  1  given  the  values  at  time  n  using  the 
following  four-step  process: 


where 


U°  =Un 

Uk  =Un  -  k  =  1..3 

V 

jjn+l  =jjA 


Oik 


1 

5-k 


V  =cell  volume 


(3.52) 


(3.53) 


The  “CFD-type”  designator  refers  to  the  fact  that  intermediate  values  of  Uk  are  over¬ 
written  each  time  step  instead  of  being  stored  at  each  step  as  in  a  traditional  Runge-Kutta 
scheme. 
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IV.  Results  of  Solutions  to  Model  Problems 


In  order  to  better  understand  the  performance  of  the  fourth-order  cRoe  scheme,  several 
model  problems  will  be  investigated.  The  first  step  is  to  verify  that  the  methodology  behind 
the  compact-difference  formulation  is  correct,  and  that  there  are  well-defined  advantages  to 
the  method.  The  next  step  is  to  duplicate  the  results  of  the  Compact  Adaptive  Interpolation 
(cRoe)  Scheme  of  Deng  and  Maekawa  |3  with  respect  to  the  one- dimensional  Euler  equa¬ 
tions.  The  cRoe  scheme  will  then  be  extended  to  two  dimensions,  and  applied  to  additional 
problems  to  understand  its  behavior  and  formal  order  of  accuracy.  Throughout  this  process 
the  cRoe  results  will  be  compared  with  the  two  methods  it  is  derived  from:  the  third-order 
Roe  (oRoe)  scheme  and  the  fourth-order  compact  (cC)  scheme. 


4.1  Linear  Wave  Equation 

As  noted  previously,  the  linear  wave  equation  (LWE)  is  a  one-dinrensional  scalar  equa¬ 
tion  that  can  be  used  to  model  the  transportive  behavior  of  more  complex  equations.  It  is 
often  a  good  starting  place  for  the  verification  of  a  new  method  due  to  it’s  simplicity.  If 
the  solution  method  works  on  the  LWE,  there  is  a  chance  that  it  will  also  work  on  a  more 
complicated  set  of  equations.  However,  if  the  method  fails  on  the  linear  problem  it  will 
mostly  likely  not  work  on  a  nonlinear  problem. 


4-1.1  Setup  Conditions.  There  are  various  explicit  solutions  to  the  LWE,  all  with 
different  strengths  and  weaknesses.  Tannehill  et  al.  [2l]  present  many  of  these  methods  along 
with  an  in-depth  analysis  of  the  solutions.  For  simplicity,  the  solution  method  considered 


here  was  the  first-order  upwind  solution  (Equation  3.7).  First,  a  one-dinrensional  domain 
of  A x/L  =  0.00625  was  used,  where  Ax  is  the  step  size  and  L  is  the  length  of  the  domain. 
The  domain  was  initialized  with  a  triangular  wave  using  the  following  initial  conditions: 


fo(x) 


0  y  <  0.075 

jL/ 

8. Ox  —  0.6  0.075  <  y  <  0.2 

—8. Ox  +  2.6  0.2  <  y  <0.325 

1j 

0  y  >  0.325 

L 


(4.1) 
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The  Courant  number  v  =  aAt/Ax  was  chosen  to  be  0.5  to  ensure  solution  stability  and 
to  demonstrate  the  dissipative  nature  of  the  first-order  upwind  scheme.  The  amount  of 
dissipation  decreases  as  v  approaches  one.  The  solution  was  allowed  to  progress  for  re  =  64 
time  steps,  which  ensured  the  traveling  wave  remained  inside  the  domain  boundaries.  The 
exact  solution  (Equation  |2.2[)  results  in  a  shift  of  the  initial  wave  a  distance  32 A x/L  to  the 
right,  with  no  decrease  in  signal  amplitude. 

This  first-order  upwind  method  was  compared  to  a  fourth-order  compact  differencing 
scheme  (Equation  |3.9[)  combined  with  a  fourth-order  Runge-Kutta  time  integration  scheme 
(Equation  |3.10[)  using  the  same  step  conditions. 


4.1.2  Linear  Wave  Equation  Results.  Figure  4.1  compares  the  results  of  the  first- 
order  and  fourth-order  methods.  Both  solutions  transported  the  wave  to  near  the  exact 
location  with  a  small  lagging  phase  error,  but  it  is  clear  that  the  compact  difference  method 
preserved  the  shape  of  the  initial  triangular  wave  much  better  than  the  first  upwind  method 
did.  The  first  upwind  method  experienced  nearly  a  30%  loss  of  maximum  amplitude  while 
the  compact  method  only  lost  5%. 


Figure  4.1:  Solutions  to  the  Linear  Wave  Equation  after  64  timesteps,  v  =  0.5 
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Figure  4.2  shows  the  error  associated  with  the  two  approximations.  Again,  it  is  clear 
that  the  first  order  wave  introduced  a  large  amount  of  dissipation  error.  The  greatest 
amount  of  error  for  both  solutions  was  near  the  peak  of  the  wave  where  the  gradients  were 
the  largest,  but  the  compact  error  was  nearly  an  order  of  magnitude  less  than  the  first 
upwind  method.  The  higher  order  compact  method  did  not  broaden  the  base  of  the  wave 
nor  substantially  decrease  the  magnitude  of  the  wave,  but  it  did  introduce  some  dispersion 
errors  (oscillations)  before  and  after  the  wave.  Both  of  these  results  can  be  traced  back 
to  the  truncation  error  terms.  The  first  order  upwind  methods  truncation  error  has  a 
dominant  even-ordered  derivative  which  causes  dissipation  effects  while  the  leading  term 
in  the  truncation  error  of  the  compact  method  is  an  odd-ordered  derivative  and  hence 
dispersive  in  nature. 


Overall,  the  fourth  order  compact  formulation  produced  a  better  approximation  of  the 
exact  solution.  The  next  step  is  to  extend  this  compact  method  to  a  nonlinear  vectorized 
set  of  equations  and  combine  it  with  Roe’s  approximate  Riemann  solver. 


Figure  4.2:  Error  in  Solutions  to  the  Linear  Wave  Equation  after  64  timesteps,  v  =  0.5 
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4-2  ID  Shock  Tube 


4-2.1  Problem  Description.  The  first  nonlinear  model  problem  considered  was  the 
one-dimensional  shock  tube  problem  introduced  by  Sod  (20  .  In  reality,  a  shock  tube  is  a 
three-dimensional  problem,  but  due  to  planar  symmetry  can  be  approximated  with  a  single 
dimension.  The  problem  consists  of  a  one-dimensional  domain  of  length  L.  The  domain 
is  split  into  two  portions,  each  of  which  have  different  initial  conditions  (noted  as  the  left 


and  right  states)  as  shown  in  Figure  4.3  The  two  portions  are  isolated  until  the  time  t  =  0 
when  they  are  allowed  to  interact.  This  problem  could  be  experimentally  considered  using 
a  gas-filled  tube  with  a  thin  membrane  separating  two  states.  The  membrane  would  be 
ruptured  at  t  =  0  and  the  two  sides  would  interact. 


Initial  Position 
Of  Separation 


i 


Figure  4.3:  Shock  Tube  Flow  State  at  t  =  0 

The  Riemann  solution  to  this  problem  allows  for  the  existence  of  a  weak  solution, 
where  the  value  and  derivatives  of  the  function  can  be  discontinuous  but  bounded  |TI] .  For 
the  shock  tube  problem,  there  exists  two  discontinuities:  a  traveling  shock  wave  due  to 
compression  and  a  contact  surface  due  to  the  different  initial  states.  The  shock  is  manifest 
in  the  solution  as  a  discontinuity  in  density,  velocity,  and  pressure.  The  contact  surface 
only  appears  as  a  discontinuity  in  density  while  pressure  and  velocity  are  continuous  and 
constant  across  it. 
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At  t  =  0.  before  the  the  two  states  are  allowed  to  interact,  the  initial  velocity  is 
uniformly  zero.  At  t  >  0  the  fluid  on  the  left  (high  pressure)  side  of  the  interface  begins  to 
move  to  the  right  (low  pressure)  side  of  the  interface.  If  the  difference  in  pressure  is  great 
enough  the  compression  waves  will  coalesce  into  a  traveling  shock  which  will  travel  to  the 
right  at  the  shock  propagation  speed.  At  the  interface  between  the  two  states  there  will 
exist  a  discontinuity  in  density  due  to  the  differences  in  initial  conditions.  This  discontinuity 
will  also  travel  to  the  right  at  a  speed  less  than  the  shock  propagation  speed.  As  fluid  in 
the  the  left  side  of  shock  tube  begins  to  travel  to  the  right  the  density  in  the  tube  will 
decrease  and  an  expansion  wave  will  travel  from  the  initial  interface  to  the  left  side  of  the 
tube.  This  expansion  wave  is  not  a  discontinuous  solution,  but  could  have  steep  gradients 
at  the  leading  and  trailing  edges. 


4- 2. 2  Exact  Riemann  Solution  to  the  Sod  Problem.  The  flow  within  the  shock 


tube  at  a  time  t  >  0  is  shown  in  Figure  4.4  The  state  of  the  system  at  t  >  0  is  described 


by  the  solution  to  the  Riemann  problem  [11] .  The  flow  can  be  divided  into  five  states:  the 
original  left  (L)  and  right  ( R )  states,  the  state  between  the  shock  and  the  contact  surface 
(2),  the  state  between  the  contact  surface  and  the  expansion  fan  (3),  and  finally  within 
the  expansion  fan  (5).  Using  Riemann  invariants,  an  implicit  relation  for  the  pressure  ratio 
P  =  P2/PR  across  the  shock  can  be  developed  in  terms  of  the  left  (L)  and  right  ( R )  states: 
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Once  the  pressure  ratio  across  the  shock  is  known  the  location  of  the  shock  at  a 
given  point  in  time  can  be  determined  using  the  propagation  speed  of  the  shock  ( C )  in  the 
undisturbed  (in  this  case  right)  region: 
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Figure  4.4:  Shock  Tube  Flow  State  at  t  >  0 


The  contact  surface  is  a  discontinuity  in  density,  but  not  in  pressure  or  velocity,  hence 
P3  =  P-2  and  W3  =  U‘2 ■  This  surface  travels  behind  the  shock  at  the  constant  speed  V 
determined  by 


V  =  U2  = 


P  -  l _ or _ u 

Vi1  +  aP )  V 7(7  -  l)/2  R 


(4.5) 


The  density  behind  the  shock  is  also  a  function  of  the  pressure  ratio  P 


P  2 


1  +  aP 


a  +  P 


PR 


(4.6) 


The  values  of  the  flow  variables  through  the  expansion  fan  (region  5)  are  determined  by  the 
left  running  characteristics  with  constant  Riemann  variables  such  that 
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Assuming  that  the  flow  is  initially  at  rest  [ul  =  ur  =  0),  the  velocity  throughout  the  fan 
varies  from  zero  to  V  described  by 


2  f(x-x0  7  -  1 

for  -  (  bj— +  a f  J"  <  ( —yt  V  -  ",  -  (4.8) 


Once  the  velocity  distribution  is  known  in  region  5  the  speed  of  sound  05  and  pressure  p§ 
can  be  expressed  in  terms  of  the  position  in  the  tube,  and  time  elapsed 


<25  —  ^5 
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4-2.3  Setup  Conditions.  For  this  problem  the  initial  conditions  were  as  follows: 


P 

u 

P 

Left 

1.0 

0.0 

1.0 

Right 

0.125 

0.0 

0.1 

The  grid  spacings  were  varied  from  A x/L  =  5E- 03  to  3.125E-04.  The  non-dimensional 
time  step  was  kept  constant  at  1.0.E-04.  The  oRoe  method  used  a  MinMod  flux  limiter, 
while  the  cRoe  scheme  used  a  k  of  0.3  (this  particular  choice  of  k  will  be  discussed  later 
in  this  section).  The  solution  was  allowed  to  progress  until  t  =  0.16  to  keep  the  shock  and 
expansion  waves  within  the  computational  domain.  Boundary  conditions  were  enforced  by 
flux- fixing  the  first  and  last  four  cells  within  the  domain  to  zero. 

The  fourth-order  cRoe  solves  systems  of  implicit  equations  twice  per  Runge-Kutta 
stage:  once  for  the  left  and  right  variable  state  interpolation,  and  then  again  to  solve  for 
the  spatial  derivatives  from  Equation [R5j  Both  of  these  steps  used  the  Thomas  tri-diagonal 
algorithm  |[2l]  to  solve  the  implicit  matrices.  The  results  of  the  cRoe  scheme  were  then 
compared  to  both  the  third-order  oRoe  scheme  and  a  first-order  Roe  scheme  (fRoe).  This 
comparison  was  made  to  show  the  dissipative  effects  of  a  first  order  solution.  The  fourth- 
order  cC  compact  scheme  was  not  considered  for  this  problem  since  it  is  unconditionally 
unstable  for  a  discontinuous  (shock)  solution. 
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4-2-4  ID  Shock  Tube  Results.  Figures  4.5  and  4.6  show  the  density  and  velocity 
distributions,  respectively,  within  the  shock  tube  at  t  =  0.16.  In  the  density  distribution  the 
three  distinct  features  of  the  flow  are  readily  visible.  The  density  at  the  far  left  and  far  right 
sides  of  the  domain  remain  at  their  initial  states,  1.0  and  0.1  respectively.  The  density  then 
decreases  through  the  expansion  wave,  and  is  then  constant  until  approximatly  x/L  =  0.65 
at  the  contact  discontinuity.  The  density  again  remains  constant  until  the  shock  near 
x/L  =  0.78.  In  the  velocity  distribution  the  left  and  right  states  remain  at  the  initial  (zero) 
condition  while  the  fluid  in  the  middle  of  the  domain  is  moving.  The  velocity  distribution 
across  the  expansion  wave  is  again  nearly  linear  while  the  shock  has  a  steep  discontinuity. 
Unlike  the  density  distribution,  the  contact  surface  is  not  visible  in  the  velocity  field  since 
normal  velocity  is  constant  across  it.  All  of  these  flow  features  are  in  agreement  with  the 
exact  solution. 


Figure  4.5:  Comparison  of  Density  Distribution  in  ID  Shock  Tube  Problem  at  t  =  0.16. 

As  expected,  the  first-order  Roe  scheme  performed  the  worst  with  a  great  deal  of 
dissipation  at  every  location  a  gradient  exists.  The  corners  near  the  expansion  wave  are 
diminished,  more  so  at  the  left  edge  than  the  right.  The  contact  surface  is  flatted  at  both 
ends  and  the  shock  is  contained  (captured)  within  17  cells  for  the  A  x/L  =  0.025  case  (see 
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Figure  4.6:  Comparison  of  Velocity  Distribution  in  ID  Shock  Tube  Problem  at  t  =  0.16. 


Figure  4.7).  Similar  to  the  behavior  noted  with  the  linear  wave  equation,  the  first  order  fRoe 


scheme  has  a  truncation  error  dominated  by  the  leading  even  derivative  terms  which  are 
dissipative.  This  combined  with  the  added  numerical  dissipation  due  to  the  Roe-averaged 
variables  results  in  a  less  than  desirable  solution. 

The  next-best  solution  was  produced  by  the  third-order  oRoe  scheme.  All  gradients  of 
density  were  captured  more  sharply  than  the  first-order  scheme,  especially  at  the  left  side  of 


the  expansion  wave  as  shown  in  Figure  4.6  Velocity  gradients  were  also  better  resolved  at 
both  the  left  and  right  sides  as  shown  in  Figure  [476}  Looking  in  detail  at  the  shock  capture 


in  Figure  4.7  the  oRoe  scheme  reduced  the  width  to  9  points  for  the  A x/L  =  0.025  case. 


All  of  these  improvements  came  from  the  fact  that  the  scheme  is  using  the  MUSCL  variable 
extrapolation  scheme  to  improve  the  overall  accuracy.  While  a  limiter  was  used,  and  hence 
the  local  order  of  accuracy  is  only  first  order  at  the  discontinuities,  the  overall  result  is  a 
better  match  to  the  exact  solution  than  the  fRoe  case. 

Finally,  the  best  solution  was  produced  by  the  compact-upwind  cRoe  scheme.  Again, 


improvements  at  all  gradients  were  made.  As  shown  in  Figures  4.5  and  4.6  the  expansion 
wave  with  the  cRoe  solution  was  very  close  to  the  Riemann  solution,  while  the  contact 
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x/L 


Figure  4.7:  Detail  of  Shock  Capture  of  Roe  Method,  Density 

surface  was  sharp.  The  shock  capture,  shown  in  Figure  |4?7|  was  further  reduced  to  5  points 
for  the  A  x/L  =  0.025  case.  Once  again  the  improvements  are  most  likely  due  to  the  higher- 
order  variable  interpolation  and  the  use  of  the  compact  differencing.  It  is  important  to  note 
that  the  cC  compact  scheme  could  not  be  used  on  this  problem  because  of  stability  issues, 
but  the  cRoe  scheme  worked  well  and  showed  no  stability  problems.  The  only  difference 
between  these  two  schemes  is  that  the  cC  scheme  uses  exact  fluxes  while  the  cRoe  scheme 
uses  Roe- averaged  fluxes.  The  results  shows  that  the  Roe  fluxes  add  enough  numerical 
dissipation  to  eliminate  the  non-physical  behavior  of  the  compact  difference  method  and 
produce  a  stable  solution. 


4-2.5  Accuracy  of  the  cRoe  Solution  Applied  to  the  ID  Shock  Tube.  To  verify  that 


Equation  3.43  holds  true  (that  the  order  of  a  compact  formulation  using  Roe  fluxes  will 
have  the  same  order  accuracy  as  the  variables  at  the  left  and  right  states)  two  methods  were 
compared:  the  third-order  oRoe  scheme,  and  a  “Compact  oRoe”  scheme  using  Equation 


1.5  with  Roe-averaged  fluxes  on  the  right-hand  side.  Figure  4.8  shows  the  overall  density 
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distribution  and  the  shock  detail  for  these  two  methods.  It  is  clear  that  there  is  little  to  no 


variation  in  the  results  of  the  two  schemes.  Additionally,  the  same  comparison  was  made 


with  the  first-order  Roe  scheme  and  eqaution  1.5  with  similar  results.  This  result  verifies 

will  only  produce  dF/dx  terms  of  same  or 


that  the  compact  formulation  of  Equation 


1.5 


lesser  order  than  the  order  of  F  that  is  used  in  the  right-hand  side  and  that  a  fourth-order 
variable  interpolation  is  needed  to  create  a  global  fourth-order  scheme. 


x/L 


Figure  4.8:  Detail  of  Normal  Shock  in  Comparison  of  Solutions  to  ID  Shock  Tube  Problem 
at  t  =  0.16. 


4-2.6  Variation  of  the  cRoe  n  Parameter  in  the  ID  shock  tube.  There  was  a  large 
variation  in  the  results  of  the  cRoe  scheme  when  the  k  switching  parameter  was  adjusted,  k 
is  a  weighting  factor  which  allows  the  compact  interpolant  to  change  bias.  In  effect  it  tells 
the  method  when  the  solution  is  non-smooth  and  when  the  stencil  should  be  biased  one 
direction  or  another.  When  k  is  too  high  the  stencil  switches  too  easily  and  unnecessary 
upwinding  can  occur.  When  k  is  zero  the  solution  reverts  back  to  a  central-difference  scheme 
which  is  unstable  across  discontinuities.  Variation  of  k  away  from  the  optimal  point  resulted 
in  oscillations  in  both  the  smooth  and  non-smooth  regions  of  the  flow.  Deng  |3  suggested 
that  k  =  0.5  was  a  good  all-around  value,  but  the  best  results  were  seen  with  n  =  0.3.  The 
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results  for  the  density  profile  around  the  contact  surface  with  variations  in  n  can  be  seen  in 


Figure  4.9  There  was  a  small  amount  of  overshoot  (oscillation)  seen  in  the  velocity  across 
the  shock  wave  which  could  not  be  removed  for  any  value  of  k. 
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Figure  4.9:  Variation  in  k  Parameter  for  ID  Shock  Tube  Problem,  Detail  of  Contact 

Surface  Discontinuity 

Deng  1 3  also  suggests  that  these  oscillations  are  caused  by  the  collision  of  the  charac¬ 
teristics  across  the  discontinuities  and  that  they  can  be  eliminated  by  performing  the  variable 
interpolations  on  the  characteristic  variables.  The  characteristic  variables  are  one  level  be¬ 
low  the  primitive  variables  and  are  the  information  that  are  carried  along  the  characteristics 
in  a  hyperbolic  solution.  Solving  on  the  characteristic  level  would  be  very  computationally 
expensive  and  was  deemed  unnecessary  for  the  problems  considered  in  this  work. 
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4-2.7  Consistency  of  the  ID  cRoe  Scheme.  As  previously  stated,  any  numerical 
discretization  of  a  PDE  is  an  approximation  and  has  an  associated  truncation  error  term. 
The  scheme  is  said  to  be  consistent  if  the  error  between  the  PDE  and  the  approximation 
vanishes  as  the  mesh  is  refined,  i.e.  Ax  — >  0.  Refining  the  mesh  to  this  point  should  make 
the  truncation  error  vanish.  If  the  scheme  is  not  consistent  it  is  not  a  useful  scheme,  as  it 
is  not  an  accurate  approximation  to  the  governing  equation.  The  general  compact  scheme 
can  be  readily  proven  to  be  consistent  through  the  substitution  of  TSEs.  The  compact 
Roe  scheme  is  more  difficult  to  prove  due  to  the  complexity  of  solving  in  terms  of  the  Roe 
averaged  variables,  and  there  consistency  is  examined  numerically. 

Figure|4.10fc,  shows  the  results  of  the  cRoe  solution  near  the  shock  as  the  grid  is  refined 
from  A x/L  =  5.0E-03  to  x/L  =  3.125E-04.  As  the  mesh  is  refined  the  shock  is  captured  in 
fewer  and  fewer  cells  and  the  steepness  of  the  shock  increases  towards  the  theoretical  value. 
It  is  apparent  from  these  results  that  the  scheme  is  consistent  as  the  solution  approaches 
the  exact  solution  rapidly  as  the  mesh  is  refined.  Since  the  scheme  appears  consistent  it  is 
a  valid  approximation  to  the  governing  equations  and  should  be  further  investigated. 


Figure  4.10:  cRoe  Solution  Consistancy  to  ID  Shock  Tube  Problem 


4-13 


4-3  2D  Shock  Tube 


Up  to  this  point,  all  of  the  results  have  been  previously  completed  by  other  parties, 
and  have  been  duplicated  here  to  ensure  repeatability  of  the  solution  methods.  The  next 
step  in  investigating  the  fourth-order  cRoe  scheme  was  to  extend  it  to  two-dimensions  and 
repeat  the  shock  tube  problem. 


4-3.1  Problem  Description  and  Setup  Conditions.  In  two-dimensions,  the  shock 
tube  problem  is  very  similar  to  the  one-dinrensional  case.  Since  the  initial  conditions  are 
only  a  function  of  one-dimension ,  the  flow  will  be  uniform  in  the  second  dimension  and  the 
results  will  be  the  same  as  the  one-dimensional  case. 


Two-dimensional  domains  with  A x/L  varied  from  5.0E-03  to  1.25.E,-03  were  con¬ 
structed.  The  domains  were  of  a  length  of  L  and  a  height  of  lOAx,  and  initialized  with  the 
same  conditions  along  the  length  as  shown  in  Table  4.2.3|  and  constant  along  the  height. 
The  same  time  steps  and  total  solution  time  as  before  were  used  for  the  two-dimensional 


case. 


4-3.2  2D  Shock  Tube  Results.  Results  for  the  overall  density  distribution  for  the 
two-dimensional  shock  tube  case  can  be  seen  in  Figure  4.11[  Comparing  this  plot  with  the 
one-dinrensional  results  in  Figure  [4~5j  it  can  be  seen  that  the  two-dimensional  results  agree 
well.  All  of  the  same  relevant  flow  features  are  present. 


4-3.3  Consistency  of  the  2D  cRoe  scheme.  Again  the  consistency  of  the  scheme  was 
verified  by  refining  the  mesh,  and  once  again  the  numerical  solution  approached  the  exact 


solution,  as  shown  in  Figure  4.12  It  was  expected  that  the  solution  would  be  consistent  in 


two  dimensions,  since  it  was  consistent  in  one  dimension,  but  this  result  verified  that  the 
method  was  properly  implemented  into  the  two-dimensional  solver. 


4-3-4  Variation  of  the  cRoe  k  Parameter  in  the  2D  Shock  Tube.  Now  that  the 
cRoe  solution  was  shown  to  be  a  consistent  solution  in  two  dimensions,  the  influence  of  the 
k  switching  parameter  was  re-checked.  Figure  4.13|r  shows  the  results  as  k,  is  varied  from 
0.1  to  0.7.  When  k  =  0  the  scheme  reverts  to  a  fourth  order  central  scheme  and  becomes 
unstable.  Genearally  speaking,  the  lower  the  value  of  n  the  better  the  discontinuities  were 
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Figure  4.11:  Comparison  of  Density  Distribution  in  2D  Shock  Tube  Problem  at  t  =  0.16. 

captured,  while  the  higher  the  value  of  n  the  more  oscillations  were  seen  in  smooth  regions 
of  the  flow.  For  very  low  values  of  n  (k  <  0.1)  there  were  some  oscillations  seen  near  the 
discontinuities.  The  results  are  similar  to  those  seen  in  the  one-dimensional  results,  most 
likely  due  to  the  fact  as  k  approaches  one  the  solution  is  more  prone  to  switching  the  bias 
of  the  interpolant  stencil.  The  more  it  switches  the  more  likely  it  will  upwind  incorrectly 
and  introduce  non-physical  behavior.  It  is  interesting  to  note  that  even  when  there  is  error 
introduced  due  to  this  switching,  the  Roe  scheme  introduces  enough  numerical  dissipation 
to  maintain  stability.  A  value  of  k  =  0.3  was  chosen  to  best  balance  the  capture  of  the 
discontinuities  and  maintain  non-oscillatory  behavior  in  the  smooth  regions. 
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Figure  4.12: 


cRoe  Solution  Consistancy  to  2D  Shock  Tube  Problem 


x/L 


Figure  4.13:  Variation  in  k  Parameter  for  2D  Shock  Tube  Problem,  Detail  of  Contact 

Surface  Discontinuity 
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4-4  Advection  of  a  Vorticial  Structure  on  a  Rectilinear  Mesh 

The  accuracy  of  a  numerical  approximation  is  determined  by  the  order  of  the  leading 
truncation  error  term.  For  a  scheme  that  is  order  n  in  space,  as  the  grid  is  refined  the 
truncation  error  should  diminish  as  Axn.  For  example,  the  theoretical  error  of  a  second 
order  solution  should  decrease  by  one- fourth  as  the  grid  is  refined  with  twice  as  many  points. 
If  the  variation  of  grid  spacing  is  plotted  against  the  solution  error  on  a  log-log  plot,  the 
resulting  line  should  have  a  slope  equal  to  the  order  of  accuracy.  This  measured  order  of 
accuracy  should  be  near  the  formal  order  of  the  scheme,  i.e.  near  the  order  of  the  leading 
truncation  error  term,  but  could  vary  due  to  other  higher-order  truncation  error  term  effects. 

To  measure  the  order  of  accuracy  of  the  cRoe  scheme  a  convecting  vortex  in  an  oth¬ 
erwise  uniform  flow  was  investigated.  This  problem  measures  the  ability  of  a  numerical 
simulation  to  accurately  advect  vortical  structures  similar  to  those  found  in  direct  and  large- 
eddy  simulations  (LES/DNS)  jl5].  Since  higher-order  methods  are  often  used  in  LES/DNS 
simulations  the  usefulness  of  the  cRoe  method  can  be  shown.  Additionally,  the  problem 
has  a  known  exact  solution  at  any  point  in  time  with  which  the  numerical  results  can  be 
compared  and  a  formal  order  of  accuracy  determination  can  be  made. 


44.1  Problem  Description.  The  vortex  chosen  has  a  central  core  of  vorticity  with 
the  same  sign  as  the  vortex  strength  surrounded  by  a  region  of  vorticity  of  the  opposite 
sign.  This  results  in  a  vortical  structure  which  has  a  total  circulation  of  zero.  The  initial 
condition  was  created  by  imposing  a  vortex  of  strength  C  at  a  location  (xc,  yc)  and  satisfying 
the  following  relations  [24]: 


U  =  Uoo  — 
Poc~P  = 


c(Vc)  exP(^r)  » 

exp(-r 2)  , 


=  exp(=f) 

r2  _  P-xA2+{y-yc)2 

1  ~  R2 


where  u,v,p  and  R  denote  the  Cartisian  velocity  components,  static  pressure,  and  vortex 
core  radius.  The  vortex  has  maximum  vorticity  at  the  center  which  decreases  to  a  minimum 
at  twice  the  core  radius.  The  vorticity  then  decreases  in  magnitude  to  zero  as  x  increases. 
The  pressure  distribution  was  obtained  by  integrating  dp/dr=pug/r  about  the  vortex  center 

0- 
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4-4-%  Exact  Solution  to  the  Vortex  Problem.  In  the  absence  of  viscous  dissipation 
the  overall  original  vortex  structure  should  remain  unchanged  in  time  transported  a  distance 
x  =  UooAt.  All  grids  and  number  of  iterations  were  chosen  such  that  the  vortex  traveled 
five  complete  trips  through  the  domain,  ending  up  in  the  same  place  that  it  started.  It  is 
important  to  note  that  the  error  in  the  solutions  have  two  major  portions:  a  phase  shift  due 
to  the  slowing  of  the  vortex  and  a  distortion  due  to  dissipation  and  dispersion.  The  phase 
shift  results  in  the  vortex  not  traveling  quite  as  far  as  it  should  have  in  theory  through  the 
domain,  while  the  distortion  is  due  to  errors  in  the  solution  method  itself.  To  be  consistent  all 
solution  results  are  measured  against  the  same  theoretical  velocity  distributions  regardless 
of  where  the  actual  solution  ended  up. 

4-4-3  Setup  Conditions.  The  nondimensional  vortex  strength  parameter  C/^u^R) 
was  set  at  0.02  to  match  previous  work  by  Visbal  and  Gaitonde  |24j  .  Constant  density  was 
assumed  initially  as  Gaitonde  suggested  that  the  uniform  density  approximation  is  suitable 
for  most  test  cases  but  could  be  calculated  using  constant  enthalpy  for  stronger  vortices. 

The  solution  was  applied  to  uniform,  rectilinear  grids  with  various  levels  of  resolution 
from  Ax/ R=Ay  /  R=0.6  to  0.1.  This  range  of  grid  spacing  was  chosen  to  both  allow  for 
variation  in  the  solution  errors  and  to  keep  total  computational  times  reasonable. 

The  non-dimensional  time  step  was  reduced  to  Att/oo/i?=0.002  which  results  in  a 
Courant-Friedrichs-Lewy  (CFL)  number  of  approximatly  0.02  for  the  finest  mesh.  It  was 
determined  numerically  that  this  time  step  allowed  the  solution  to  be  sufficiently  time- 
resolved  at  the  finest  grid  spacings  which  would  mean  it  was  also  time  resolved  on  the 
coarsest  mesh. 

For  the  oRoe  solution  the  minmod  flux  limiter  was  removed  to  allow  a  global  third- 
order  solution.  Limiters  work  by  reverting  to  first  order  near  non-snrooth  regions.  Even 
when  the  solution  is  “mostly  smooth”  there  still  is  possibility  of  limiter  influence  which 
would  degrade  the  accuracy  of  the  solution.  Since  the  problem  is  continuous  the  use  of  a 
limiter  is  unneeded  and  it  was  subsequently  removed. 

4-4-4  Periodic  Boundary  Condition.  Since  the  vortex  is  advected  by  the  free- 
strearn  flow  the  grid  would  have  to  be  made  very  large  or  the  solution  would  have  to  be 
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Mesh  1 


Mesh  2 


Figure  4.14:  Mesh  Overlap  on  Periodic  Boundaries 

allowed  to  leave  one  portion  of  the  domain  and  re-enter  at  another.  Again,  to  reduce 
computational  times  the  latter  was  chosen  and  all  four  sides  of  the  domain  were  treated 
with  periodic  boundary  condition.  A  five  point  overlap  was  chosen  as  shown  in  Figure  [4 .14 
such  that  the  traveling  vortex  was  allowed  to  exit  one  side  of  the  domain  and  re-enter  the 
opposite  side  in  a  cyclic  manner.  Likewise,  any  velocity  perturbations  in  the  top  and  bottom 
of  the  mesh  were  allowed  to  exit  and  re-enter  the  domain.  Visbal  and  Gaitonde  [5j  have 
shown  five  points  are  a  sufficient  overlap  to  preserve  the  accuracy  of  the  solution  on  the 
periodic  boundary  and  little  degradation  of  the  solution  occurs. 

Mesh  dimensions  of  —6  <  x  <  6  and  —6  <  y  <  6  were  chosen  such  that  at  the  initial 
state  the  velocity  at  the  boundaries  of  the  domain  were  near  machine  zero.  A  smaller  domain 
could  have  allowed  feedback  across  the  domain  in  every  iteration,  whereas  the  larger  domain 
only  transmitted  information  across  the  boundary  when  the  vortex  had  traveled  across  the 
domain. 

The  implicit  calculation  of  the  convective  gradients  ( F '  and  G ')  and  the  implicit 
calculation  of  the  left  and  right  variable  states  for  the  cRoe  method  were  modified  due  to 
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the  periodic  boundaries.  The  implicit  matrix  was  now  considered  as  a  full  periodic  tri¬ 
diagonal  system  where  the  first  and  last  rows  of  the  left-hand  side  contain  a  non-zero  term 
in  the  last  and  first  column  respectively.  This  resulting  system  was  directly  solved  using  a 
periodic  tri-diagonal  solver  [10  . 


4-4-5  Results  to  the  Vortex  Problem.  Figure  4.15  shows  contours  of  constant 
vorticity  magnitude  at  t  =  8.0  for  the  exact,  oRoe,  cC,  and  cRoe  methods  on  the  A x/R  =  0.2 
mesh.  It  is  clear  that  all  methods  did  fairly  well  preserving  the  overall  velocity  distribution 
near  the  center  of  the  domain,  but  all  methods  suffered  towards  the  outer  edges  of  the 
vortex.  The  oRoe  method  showed  a  maximum  error  along  planes  7t/4  from  the  coordinate 
axes  consistent  with  published  results  |24j.  Gaitonde  noted  that  this  is  due  to  a  bias  inherit 
to  the  Roe  scheme  itself  and  is  generally  harmless  for  most  solutions  |6|.  The  fourth- 
order  cRoe  and  the  fourth-order  cC  schemes  both  show  similar  results,  with  slightly  better 
flow  preservation  than  the  third-order  oRoe  scheme.  All  three  methods  continued  to  act 
consistently  as  they  approached  the  theoretical  value  as  the  mesh  was  refined. 

Figures  |4.16  4.17  and  4.18  show  the  swirl  velocity  and  vorticity  profiles  along  the 
y  =  0  centerline  for  the  oRoe,  cC,  and  cRoe  methods  respectively.  The  oRoe  method 
showed  the  most  dissipation  error  for  all  values  of  A  x/R  with  the  largest  errors  on  the 
coarsest  mesh.  For  A  x/R  =  0.3,  the  resulting  peak  vorticity  magnitude  was  only  40%  of 
the  initial  value.  This  result  is  contrasted  to  the  cC  and  cRoe  schemes  which,  on  this  mesh, 
recovered  93%  and  86%,  respectively,  of  the  original  peak  values.  The  oRoe  scheme  also 
caused  a  smoothing  of  the  velocity  distribution  with  the  peak  values  diminished  and  the 
vortex  width  broadened,  both  due  to  the  inherit  dissipative  nature  of  the  oRoe  scheme.  It 
is  worth  noting  that  a  first-order  Roe  scheme,  similar  to  the  fRoe  shown  on  the  shock  tube, 
would  have  completely  dissipated  the  vortex  after  one  cycle  through  the  domain. 

The  cC  compact  method  did  the  best  job  of  the  three  schemes  at  preserving  the  peak 
velocity  values,  but  was  also  the  most  dispersive  of  the  three.  For  many  of  the  cC  test  cases 


shown  in  Figure  4.17  the  solution  exhibited  a  lagging  phase  error  with  the  vortex  slightly 
shifted  to  the  left  of  the  theoretical  location,  and  more  so  on  the  coarsest  mesh.  There  was 
also  a  large  amount  of  oscillations  at  the  trailing  edge  of  the  vortex,  also  due  to  dispersion. 
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Figure  4.15:  Vorticity  Profiles  of  Advecting  Vortex  at  T=30.0,  Ax/ii=0.20 
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As  mentioned  previously,  the  fourth-order  scheme  has  a  truncation  error  dominated  by  an 
odd-ordered  derivative  term  which  produces  the  dispersion  effect. 


The  cRoe  method,  (Figure  4.18),  exhibited  similar  results  to  both  the  oRoe  and  cC 
schemes  (Figures  4.16  and  Em  respectively),  most  likely  since  it  contains  elements  of  each 
of  them.  The  cRoe  scheme  did  a  better  job  of  preserving  peak  values  of  velocity  than  the 
oRoe  method,  but  not  quite  as  good  a  job  as  the  cC  compact  scheme  did.  However,  the 
cRoe  scheme  did  not  display  as  much  dispersion  or  phase  error  as  the  cC  scheme  did,  most 
likely  due  to  the  added  numerical  dissipation  from  the  Roe  fluxes. 

Figures  |4.19  4.20,  and  4.21  show  the  associated  error  of  the  oRoe,  cC,  and  cRoe 
methods  respectively,  for  both  swirl  velocity  and  vorticity  magnitude  along  the  y  =  0 
centerline.  Again,  the  oRoe  method  had  the  largest  errors  associated  for  the  velocity  profile, 
with  errors  of  nearly  an  order  of  magnitude  more  than  the  cC  and  cRoe  compact  schemes. 
Overall  the  cC  scheme  had  the  lowest  overall  absolute  errors  based  on  a  L ^  norm,  with  the 
cRoe  coming  in  second,  and  the  oRoe  last. 


The  oRoe  errors  were  smoothly  distributed  across  the  domain  with  a  shape  similar 
to  the  associated  function.  The  compact  schemes  produced  errors  of  a  lesser  magnitude, 
but  with  more  variation  across  the  domain,  similar  to  the  errors  seen  in  the  linear  wave 
equation.  Again,  this  could  be  attributed  to  the  numerical  dissipation  of  the  oRoe  scheme 
and  the  dispersive  nature  of  the  compact  difference.  The  compact  schemes  do  a  better 
job  of  resolving  the  high-frequency  waves  in  the  solution,  whereas  the  oRoe  scheme  damps 
them  out.  This  makes  the  oRoe  scheme  slightly  more  stable  and  robust  than  the  compact 
schemes,  but  also  removes  small  scale  information  from  the  flow.  If  the  goal  of  the  solution 
is  to  produce  a  well-behaved,  stable  result  the  oRoe  scheme  would  be  a  good  choice.  If 
the  goal  of  the  solution  is  to  capture  as  much  of  the  flow  structure  as  possible,  with  the 
possibility  of  some  instabilities,  the  compact  methods  would  be  best. 
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Figure  4.16: 
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Figure  4.17:  Comparison  of  cC  Solutions  to  Vortex  Problem  along  centerline  after  t  =  30. 


Figure  4.18: 
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Comparison  of  cRoe  Solutions  to  Vortex  Problem  along  centerline  after 
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Figure  4.20:  Comparison  of  Error  in  cC  Solutions  to  Vortex  Problem  along  centerline 

t  =  30. 
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4-4-6  Order  of  Accuracy  Determination.  For  the  order  of  accuracy  determination, 
an  L2  norm  of  the  swirl  velocity  along  the  y  =  0  centerline  was  calculated  and  plotted 
against  the  non-dimensional  grid  spacing  Ax / R.  The  swirl  velocity  was  chosen  because  it  is 
not  a  derived  quantity  (like  vorticity),  the  use  of  which  could  inadvertently  add  additional 
error  to  the  solution.  The  L2  norm  was  chosen  as  it  provided  the  best  match  to  published 
results  for  the  third-order  oRoe  scheme  and  the  fourth-order  cC  scheme.  The  computed 
orders  of  accuracy  based  on  a  least-squares  fit  of  the  log  of  the  error  versus  the  log  of  the 


grid  spacing  are  presented  graphically  in  Figure  4.22  and  numerically  Table  4.1 


Table  4.1:  Theoretical  and  Calculated  Orders  of  Accuracy  for  Various  Solution  Methods 
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Figure  4.22:  Order  of  Accuracy  of  Solution  Methods 
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The  orders  of  accuracy  of  the  oRoe  method  and  the  cC  method  were  comparable  with 
both  the  theoretical  and  published  values  [24].  The  cRoe  method  was  slightly  lower  than 
anticipated,  as  it  was  expected  to  be  a  fourth-order  solution.  However,  when  the  distinction 
between  order  of  accuracy  and  absolute  error  is  made,  the  cRoe  method  clearly  has  a  large 
advantage  over  the  oRoe  scheme.  For  a  consistent  scheme,  order  of  accuracy  says  nothing 
about  how  close  a  method  is  to  approximating  an  exact  value,  but  rather  the  rate  at  which 
the  solution  approaches  the  exact  value.  For  the  case  of  the  cRoe  scheme,  while  it  is  not 
approaching  the  exact  solution  at  as  great  a  rate  as  the  cC  compact  scheme,  it  produced 
absolute  errors  comparable  to  those  of  the  cC  scheme  and  much  lower  than  the  oRoe  scheme. 
The  fact  that  the  error  was  comparable  with  the  cC  method,  but  that  it  could  be  applied 
to  a  discontinuous  solution  is  a  very  important  finding. 

While  the  method  presented  for  determining  order  of  accuracy  is  theoretically  sound, 
it  does  have  it’s  pitfalls.  For  highly  accurate  schemes  there  is  a  delicate  balance  between 
producing  too  good  of  a  solution  with  a  highly  machined  mesh  and  producing  a  very  poor 
solution  due  to  a  coarse  mesh.  The  oRoe  method  was  run  on  more  refined  meshes  than  the 
compact  schemes,  principally  due  to  this  fact.  When  the  compact  schemes  were  run  with  a 
mesh  spacing  less  than  A x/R  =  0.2  the  error  calculations  decreased  in  slope  drastically,  or 
leveled-off  altogether.  Physically  what  is  happening  is  that  at  these  levels  of  mesh  refinement 
the  compact  methods  preserve  the  flow  so  well  that  the  errors  have  reached  the  range  of 
machine  zero  and  cannot  be  improved  upon.  The  only  option  for  this  case  is  to  coarsen  the 
mesh  and  try  to  increase  the  error  of  the  solution.  However,  when  the  mesh  is  coarsened 
the  oRoe  method  had  the  opposite  problem.  The  oRoe  method  was  unable  to  produce  a 
solution  on  a  mesh  with  spacing  greater  than  A  x/R  =  0.4  due  to  the  magnitude  of  the  error 
terms.  When  the  mesh  is  very  coarse  the  solution  error  behaves  very  nonlinear ly  and  it  is 
difficult  to  draw  conclusions.  The  grid  with  A  x/R  =  0.3  proved  to  be  the  point  in  common 
which  all  solution  methods  worked  well.  It  is  this  problem  of  balancing  the  mesh  refinement 
that  causes  the  orders  of  accuracy  to  be  slightly  different  than  the  theoretical  values.  These 
problems  are  further  aggravated  by  the  fact  that  the  advecting  vortex  is  a  relatively  simple 
problem.  In  order  to  introduce  more  error  for  the  compact  methods,  a  more  difficult  problem 
would  be  needed.  The  drawback  to  attempting  a  more  difficult  problem  is  that  they  often 
have  less  defined  exact  solutions  for  comparison.  Again,  even  though  the  three  methods 
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only  shared  one  common  mesh  refinement,  it  is  clear  that  the  compact  methods  produced 
much  lower  absolute  error  than  the  oRoe  method  did.  This  result  means  that  a  solution 
could  be  run  on  a  coarser  mesh  with  the  compact  methods  than  the  oRoe  method  but  still 
have  a  lower  overall  absolute  error. 
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4-5  Advection  of  a  Vortical  Structure  on  a  Curvilinear  Mesh 

Thus  far  two  dimensional  results  have  been  presented  on  uniform  rectilinear  meshes. 
While  these  meshes  are  suitable  for  analytic  investigations,  they  are  not  realistic  problems. 
There  are  not  too  many  “real-world”  problems  that  occur  with  flow  on  a  uniform  rectilinear 
domain.  In  order  to  be  useful,  a  solution  must  be  able  to  handle  variations  in  the  mesh  due  to 
stretching  and  skewing.  On  a  two-dimensional  rectilinear  mesh  the  physical  domain  exactly 
aligns  with  the  computational  domain  and  the  associated  grid  metrics  are  either  unity  or 
zero.  When  the  metrics  are  non-zero  they  can  have  a  great  influence  on  the  outcome  (and 
accuracy)  of  the  solution. 


4-5.1  Problem  Description  and  Setup  Conditions.  To  understand  how  the  solution 
methods  perform  on  a  smooth  but  non-orthogoal  curvilinear  grid  a  new  set  of  meshes  were 
generated  analytically  according  to  |5] : 


Xi,j  =xmin  +  A x0[(i  -  1)  +  Axsin(nxn(j  -  1)A y0/Ly  +  j(fx/ii) 
Vi,j  =ymin  +  A yQ[(j  -  1)  +  Aysin(nyn(i  -  1)A x0/Lx  +  jcj>y/jj ) 
Ax0  =Lx/{ii  +  1) 

A  y0  =Ly/(jj  +  1) 

1  <  i  <  ii 
1  <  j  <  33 


(4.10) 


ii  and  jj  denote  the  number  of  grid  points  in  the  x  and  y  directions,  respectively,  with 
Lx  =  xmax  -  xmin  and  Ly  =  ymax  -  ymin.  The  amplitude,  frequency,  and  phase  shift 
parameters  were  specified  as  Ax  =  Ay  =  0.5,  nx  =  ny  =  4,  and  4>x  =  4>y  =  0.0,  respectively. 
The  lengths  and  number  of  nodes  were  set  such  that  the  average  node  spacing  matched  the 
previous  rectilinear  results.  The  resulting  mesh  has  a  maximum  aspect  ratio  of  1.8,  a  ratio 
of  maximum  to  minimum  cell  size  of  1.84,  and  a  maximum  deviation  from  orthogonality 


(equi-angle  skewness)  of  0.42.  An  example  mesh  is  shown  in  Figure  4.23 


The  three  solution  methods  were  used  to  solve  the  same  advecting  vortex  problem 


described  by  equation  4.10  with  the  same  solution  settings  and  methods.  The  grid  metrics 
which  map  the  curvilinear  domain  to  the  computational  domain  were  computed  using  a 
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5-point  Overlap 


Figure  4.23:  A x/R  =  0.2  Curvilinear  Mesh  with  5-point  Overlap 


fourth-order  compact  difference  method  similar  to  Equation  3.3  to  ensure  a  global  higher- 
order  accurate  solution. 


4-5.2  Curvilinear  Vortex  Advection  Results.  Figures  4.24,  4.25  and  4.26  show  the 
centerline  swirl  velcocity  and  vorticity  magnitudes  along  y  =  0  for  the  oRoe,  cC,  and  cRoe 
methods,  respectively.  Again,  as  somewhat  expected,  the  oRoe  scheme  did  the  worst  with 
a  great  deal  of  dissipation  in  all  test  cases.  After  five  trips  through  the  domain  the  oRoe 
method  had  retained  less  than  40%  of  the  initial  signal.  The  cC  and  cRoe  methods  did 
much  better  than  the  oRoe  method,  but  still  had  greater  dissipation  and  dispersion  errors 
than  the  rectilinear  results.  This  was  an  anticipated  result,  as  the  curvilinear  mesh  is  a 
more  difficult  problem  to  solve  and  should  introduce  more  error  into  the  solution.  As  seen 
in  other  problems,  this  error  is  manifest  mainly  as  dissipation  in  the  oRoe  method  and  as 
dispersion  in  the  compact  methods. 


When  the  cC  and  cRoe  methods  are  compared  on  Figures  4.25  and  4.26,  respectively, 
the  cRoe  appeared  to  produce  better  results  on  the  curvilinear  mesh  than  the  straight 
compact  method  did.  The  cC  method  had  a  large  amount  of  oscillations  on  the  middle 


4-33 


A x/R  =  0.3  mesh,  while  the  cRoe  scheme  was  relatively  nonoscillatory.  This  could  be 
attributed  to  the  added  numerical  dissipation  the  Roe  scheme  adds  to  the  solution,  removing 
the  oscillations. 


Figures  4.27  4.28  and  4.29  show  the  error  distribution  for  swirl  velocity  and  vorticity 
magnitude  for  the  oRoe,  cC,  and  cRoe  methods  respectively  along  the  y  =  0  centerline. 
Results  were  similar  to  those  of  the  rectilinear  mesh  for  all  schemes,  except  the  magnitudes 
of  error  nearly  doubled.  While  the  L0 0  norms  for  the  cC  solutions  were  almost  always  lower 
than  the  cRoe  solutions,  the  cRoe  solution  had  much  less  dispersion  error  and  seemed  to  do 
a  better  overall  job  of  preserving  the  initial  conditions. 
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(a)  Swirl  Velocity 


x/R 


(b)  Vorticity  Magnitude 

Figure  4.24:  Comparison  of  oRoe  Solutions  to  Vortex  Problem  on  a  Curvilinear  Mesh 

Along  Centerline  After  t  =  30. 
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(b)  Vorticity  Magnitude 

Figure  4.25:  Comparison  of  cC  Solutions  to  Vortex  Problem  on  a  Curvilinear  Mesh  Along 
Centerline  After  t  =  30. 
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(a)  Swirl  Velocity 
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(b)  Vorticity  Magnitude 

Figure  4.26:  Comparison  of  cRoe  Solutions  to  Vortex  Problem  on  a  Curvilinear  Mesh 

Along  Centerline  After  t  =  30. 
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(a)  Swirl  Velocity  Error 
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(b)  Vorticity  Magnitude  Error 


Figure  4.27:  Comparison  of  Error  in  oRoe  Solutions  to  Vortex  Problem  on  a  Curvilinear 

Mesh  Along  Centerline  After  t  =  30. 
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(a)  Swirl  Velocity  Error 
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(b)  Vorticity  Magnitude  Error 


Figure  4.29:  Comparison  of  Error  in  cRoe  Solutions  to  Vortex  Problem  on  a  Curvilinear 

Mesh  Along  Centerline  After  t  =  30. 
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4-5.3  Order  of  Accuracy  Determination  on  Curvilinear  Mesh.  Again,  the  order  of 
accuracy  was  determined  with  a  least-squares  fit  to  the  error  and  grid  spacing  data.  The 
results  of  these  calulcations  are  shown  in  table  14.21 

Table  4.2:  Theoretical  and  Calculated  Orders  of  Accuracy  for  Various  Solution  Methods 
on  a  Curvilinear  Mesh  _ 


oRoe 

cC 

cRoe 

Theoretical 

3.0 

4.0 

4.0 

Computed 

2.4 

3.3 

3.8 

Ax/R 


Figure  4.30:  Order  of  Accuracy  of  Solution  Methods  on  a  Curvilinear  Mesh 

The  order  accuracy  of  the  oRoe  and  cC  methods  decreased  slightly,  which  was  expected 
and  in  accordance  with  published  results  [24.  .  All  of  the  methods  investigated  are  derived 
using  some  sort  of  difference  equation,  either  a  finite  difference  or  a  compact  difference.  Both 
of  these  difference  methods  assume  a  constant  node  spacing  Ax  in  their  derivation  using 
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Taylor  series  expansions.  On  the  rectilinear  mesh  this  is  correct,  but  it  is  not  true  on  the 
curvilinear  mesh.  The  cell  metrics  take  some  of  this  non-uniformity  into  account,  but  they 
do  not  eliminate  its  effect  all  together.  Additionally,  the  order  of  accuracy  determination 
assumes  a  constant  node  spacing,  and  that  all  cells  change  at  the  same  rate  as  the  mesh 
is  refined.  For  a  rectilinear  mesh  this  is  also  correct,  but  in  a  mesh  with  skewed  cells  the 
spacing  is  non-uniform  and  cell  sizes  will  change  at  different  rates.  However,  the  results  of 
an  order  of  accuracy  determination  on  a  curvilinear  mesh  is  still  interesting  for  it  shows  the 
utility  of  a  solution  on  a  non-ideal  mesh.  Since  most  physical  geometries  require  non-ideal 
grids,  it  is  a  worthwhile  experiment. 

The  surprise  for  this  experiment  was  the  fact  that  the  slope  of  the  cRoe  error  curve 
increased  over  that  of  the  rectilinear  results.  The  same  cautionary  statement  applies  to  the 
curvilinear  results  that  did  to  the  rectilinear  results:  the  order  of  accuracy  is  only  a  measure 
of  the  rate  of  change  of  error,  while  the  absolute  error  is  generally  the  more  useful  result. 
The  order  of  accuracy  dictates  the  rate  at  which  the  solution  approaches  the  exact  value, 
while  the  absolute  error  is  the  difference  between  the  approximation  and  the  exact  value. 
If  the  absolute  error  is  already  low,  the  rate  at  which  it  approaches  the  exact  solution  is 
not  as  important.  For  the  curvilinear  meshes  the  cRoe  scheme  constantly  produced  results 
with  less  absolute  error  than  the  oRoe  scheme,  and  sometimes  the  cC  scheme.  The  fact 
that  the  cRoe  scheme  absolute  error  is  comparable  to  the  cC  scheme  and  it  can  be  used  on 
discontinuous  solutions  is  a  very  promising  result. 


4-5-4  Variation  of  cRoe  n  Parameter  in  Vortex  Problem.  As  seen  with  the  shock 
tube  problem,  the  k  weighting  factor  had  an  influence  on  the  performance  of  the  scheme. 
One  can  recall  that  this  k  factor  is  a  switch  to  tell  the  compact  adaptive  interpolation  which 
direction  to  go.  A  value  of  k  =  0  gives  the  interpolant  a  fixed  central  stencil,  while  n  >  0 
allows  the  interpolant  to  switch  directions.  Larger  values  of  k  mean  the  switch  is  more 
sensitive  to  the  smoothness  of  the  solution. 


Figure  4.31  shows  the  results  on  a  mesh  with  A x/R  =  0.4  at  t  =  1.0,  with  k  varied 
from  0.0  to  1.0.  As  k  was  increased  from  zero  numerical  instabilities  began  to  appear  near 
the  core  radius  of  the  vortex.  As  k  increased  these  instabilities  grew  larger  and  larger. 
It  is  interesting  to  note  that  the  oscillations  were  always  bounded  and  the  solution  never 
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K=0.0 


K=0.1  K=0.3 


Figure  4.31:  Variation  in  k  parameter  for  cRoe  problem,  A x/R  =  0.4 

diverged,  but  by  n  =  0.3  the  solution  had  degraded  quite  a  bit.  The  oscillations  are  most 
likely  due  to  the  scheme  performing  an  incorrect  upwind  interpolation  due  to  an  incorrect 
switching  function.  As  k  increases  the  switching  function  becomes  more  prone  to  shifting 
to  the  left  or  right,  even  when  the  solution  was  smooth  enough  to  remain  in  a  central  mode. 
The  only  value  of  k  that  produced  acceptable  results  was  that  of  k  =  0,  which  essentially 
turns  off  the  upwind  interpolation.  This  result  is  somewhat  expected  because  the  solution  is 
smooth  and  should  never  need  to  use  an  interpolant  other  than  a  central  stencil.  However, 
it  also  shows  that  the  switching  function  that  determines  k  is  not  well  formed,  and  needs  to 
be  refined  to  better  measure  the  smoothness  of  the  stencil.  If  this  switching  function  were 
correct,  the  scheme  would  always  use  a  central  interpolant  on  the  vortex  problem  regardless 
of  the  value  of  k  chosen. 
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V.  Conclusions 


A  first-order  first  upwind  method  and  fourth-order  compact  difference  method  were  applied 
to  the  linear  wave  equation.  The  results  not  only  demonstrated  that  the  compact  difference 
is  a  valid  discretization  method,  but  showed  that  the  accuracy  and  resolving  ability  of  the 
compact  difference  is  much  greater  than  a  first-order  method. 

The  formulation  of  the  compact  difference  was  then  modified  to  use  Roe  averaged 
fluxes  on  the  right-hand  side  and  was  applied  to  the  one-dinrensional  shock  tube  problem. 
The  results  confirmed  that  the  order  of  accuracy  of  the  results  matched  the  order  of  the  left 
and  right  variable  states  used  in  the  flux  calculation. 

The  Compact  Adaptive  Interpolation  scheme  of  Deng,  et.  al.  |3|,  or  cRoe  scheme, 
was  then  incorporated  into  the  one-dimensional  solver  and  applied  to  the  one-dimensional 
shock  tube.  Previous  results  were  duplicated,  and  it  was  shown  that  the  cRoe  method  could 
produce  sharper  shock  captures  and  closer  results  to  the  exact  solution  than  the  first-order 
fRoe  and  third-order  oRoe  schemes.  The  shock  tube  mesh  was  refined  and  it  was  shown 
that  the  cRoe  scheme  behaved  consistently,  approaching  the  exact  solution  as  Ax  decreased. 
It  was  also  shown  that  the  value  of  the  cRoe  switch  weighting  parameter  k  needed  to  be 
tuned  to  produce  the  best  solution,  k  =  0.3  appeared  to  produce  the  best  shock  captures 
with  the  least  oscillations  or  overshoots. 

The  cRoe  scheme  was  then  extended  to  two- dimensions  and  included  in  a  generalized 
coordinate  solver.  The  scheme  was  then  applied  to  a  two-dimensional  shock  tube  with 
results  similar  to  the  one-dinrensional  case.  Again,  it  was  shown  that  the  method  was 
consistent  and  approached  the  exact  solution  as  the  mesh  was  refined.  Additionally,  the  k 
switch  parameter  was  varied  and  shown  to  have  an  influence  on  the  quality  of  the  solution. 

The  cRoe  scheme  was  then  applied  to  a  smooth  problem  involving  the  advection  of  a 
vortical  structure  on  rectilinear  and  curvilinear  meshes.  Results  were  compared  to  those  of 
the  third-order  oRoe  scheme  and  the  fourth-order  cC  scheme.  The  formal  order  of  accuracy 
of  the  cRoe  scheme  was  shown  to  lie  somewhere  between  the  third  and  fourth  order,  with 
an  order  of  3.4  on  the  rectilinear  mesh  and  3.8  on  the  curvilinear  mesh.  Variation  of  the  k 
parameter  showed  that  any  value  of  k  greater  than  zero  produces  undesirable  results  for  a 
smooth  problem. 
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5.1  Strengths  of  the  cRoe  scheme 

The  first  important  conclusion  with  regards  to  the  cRoe  scheme  is  that  it  is  a  consistent 
scheme  in  one  and  two  dimensions,  and  therefore  is  a  valid  approximation  to  the  Euler 
equations.  If  the  scheme  were  not  consistent  it  would  not  be  able  to  model  the  physical 
behavior  described  by  the  governing  equations.  Although  the  formal  order  of  accuracy  of  the 
cRoe  scheme  appears  to  be  slightly  less  than  that  of  the  cC  compact  scheme,  the  absolute 
error  is  comparable.  Both  the  order  of  accuracy  and  absolute  error  of  the  cRoe  scheme  were 
better  than  that  of  the  third-order  oRoe  scheme. 

The  real  strength  of  the  cRoe  scheme  is  that  it  has  an  absolute  error  comparable  to  a 
standard  compact  difference  method,  but  it  can  be  applied  to  solutions  with  discontinuities. 
The  performance  of  the  cC  compact  scheme  was  not  measured  on  the  one  or  two  dimensional 
shock  tube  problems  because  the  method  is  unconditionally  unstable  near  discontinuities. 
The  cRoe  scheme  was  able  to  produce  better  results  than  the  oRoe  scheme,  and  do  so 
without  the  use  of  limiters  or  filters. 

5.2  Weaknesses  of  the  cRoe  scheme 

The  quality  of  the  results  from  the  cRoe  scheme  depended  greatly  on  the  value  of  the 
k  switching  parameter.  There  was  a  slight  dependence  on  the  value  of  n  when  the  solution 
had  a  discontinuity,  but  a  large  dependence  when  the  solution  was  smooth.  The  fact  that 
the  interpolation  needed  to  be  set  to  a  pure  central  difference  when  the  solution  was  smooth 
was  expected,  but  the  switching  method  of  the  scheme  should  have  been  able  to  do  this  on 
it’s  own  without  setting  n  =  0. 

While  the  cRoe  method  did  improve  the  accuracy  of  the  solution  over  the  oRoe  case, 
it  also  added  a  considerable  number  of  computations  per  time  step.  For  each  conserved 
variable,  the  cRoe  scheme  must  solve  two  implicit  systems  across  the  entire  domain:  one 
for  the  left  and  right  variable  states,  and  another  for  the  flux  derivative  calculations  each 
time  step.  These  implicit  systems  are  large  systems  of  linear  equations  which  can  be  com¬ 
putationally  expensive  to  solve.  The  cRoe  scheme  always  took  the  longest  time  to  solve  the 
advecting  vortex  problem.  The  cC  scheme  was  second  and  the  oRoe  method  the  fastest. 
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For  a  large  domain  with  many  grid  points  the  difference  in  number  of  computations  could 
increase  the  solution  time  dramatically. 

5.3  Relative  Computational  Costs 

The  cC  and  cRoe  schemes  had  a  higher  computational  cost  than  the  oRoe  scheme 
due  to  the  implicit  calculations  involved.  The  cC  scheme  had  to  solve  the  implicit  flux 
calculation  each  Runge-Kutta  stage,  while  the  cRoe  method  solved  both  the  implicit  variable 
state  problem  and  the  implicit  flux  problem  each  Runge-Kutta  stage.  However,  the  order 
of  accuracy  determinations  proved  that  the  cC  and  cRoe  methods  had  a  lower  absolute 
error  than  the  oRoe  method.  This  means  the  cC  and  cRoe  methods  could  be  applied  to 
a  coarse  grid  and  produce  errors  comparable  to  the  oRoe  method  applied  to  a  fine  grid. 
Numerically  it  was  shown  that  the  oRoe  method  on  the  Ax/R=0.2  mesh  produced  similar 
results  as  the  cC  and  cRoe  methods  on  the  Ax/R=0A  mesh.  For  a  complex  problem,  it  is 
conceivable  the  additional  expense  of  the  cC  and  cRoe  methods  could  be  offset  due  to  the 
reduced  number  of  mesh  points.  This  was  not  the  case  with  the  advecting  vortex  results,  as 
the  solution  times  for  the  cC  and  cRoe  methods  took  11%  and  18%  longer  than  the  oRoe 
cases,  respectively. 

5-4  Future  Work 

The  next  step  to  validating  this  new  cRoe  scheme  would  be  to  apply  it  to  the  full 
Navier-Stokes  equations  and  use  it  to  solve  a  viscous  problem  such  as  a  flat  plate  boundary 
layer.  The  numerical  results  could  then  be  compared  to  the  Blasius  solution  for  accuracy. 
The  addition  of  viscosity  would  add  another  level  of  difficulty  to  the  solution.  Since  viscous 
effects  occur  on  a  scale  much  smaller  than  convective  effects  the  resolution  ability  of  the 
compact  schemes  could  be  investigated.  In  order  to  maintain  a  global  order  of  accuracy  the 
viscous  fluxes  of  the  Navier-Stokes  equations  should  be  calculated  to  the  same  or  greater 
order  of  accuracy  as  the  convective  fluxes.  For  the  cC  and  cRoe  schemes  this  could  be 
accomplished  using  similar  compact  difference  formulations. 

Another  important  problem  that  should  be  considered  would  be  one  involving  a  stand¬ 
ing  discontinuity,  such  as  a  high  Mach  number  blunt-body  or  a  ramp  with  an  oblique  shock. 
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The  shock  tube  problem  demonstrated  the  ability  of  the  cRoe  method  to  handle  a  discon¬ 
tinuous  solution  on  a  transient  problem.  It  is  often  more  difficult  for  numerical  solutions  to 
resolve  a  stationary  discontinuity  than  a  moving  one  without  excessive  smearing  or  broad¬ 
ening.  For  error  comparison  purposes  there  exists  experimental  data  on  blunt  bodies  and 
an  exact  analytical  solution  for  oblique  shock  problems. 

Finally,  any  refinements  that  could  be  made  to  the  adaptive  interpolation  method 
would  be  an  improvement.  For  the  method  to  be  more  robust  the  interpolation  should  be 
able  to  better  gauge  the  smoothness  of  the  function  without  a  manual  adjustment  of  the 
value  of  k.  Several  flux  splitting  schemes,  such  as  Edward’s  Low  Diffusion  Flux  Splitting 
Scheme  [4],  have  logic  incorporated  into  them  to  determine  sonic  transitions  and  contact 
surfaces  for  proper  upwinding.  If  similar  logic  could  be  included  in  the  cRoe  method,  the 
robustness  of  the  scheme  would  be  greatly  increased. 

5.5  Final  Conclusion 

The  cRoe  scheme  proposed  by  Deng  and  Maekawa  |3|  is  a  consistent  compact  dif¬ 
ference  scheme  with  near  fourth-order  accuracy.  It  was  shown  to  capture  discontinuities 
and  gradients  in  a  shock  tube  problem  better  than  a  third-order  Roe  scheme,  and  was  also 
able  to  produce  absolute  errors  similar  to  that  of  a  standard  compact  scheme  on  a  smooth 
problem.  While  its  robustness  could  be  improved  through  a  better  interpolant  switching 
parameter,  it  successfully  produced  a  nonoscillatory  result  on  a  discontinuous  solution  using 
an  unfiltered  compact  difference  method. 
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and  Test  Engineer  for  the  Advanced  Tactical  Laser  1-rnol  Demonstrator. 

In  August  of  2003,  he  entered  the  Aeronautical  Engineering  program  at  the  Graduate 
School  of  Engineering  and  Management,  Air  Force  Institute  of  Technology.  His  academic 
focus  areas  were  Computational  Fluid  Dynamics  and  Aerodynamics.  Upon  graduation,  he 
will  be  assigned  to  the  Computational  Research  Branch,  Aeronautical  Sciences  Division,  Air 
Vehicles  Directorate,  Air  Force  Research  Lab  at  Wright-Patterson  AFB,  OH  where  he  will 
continue  his  work  in  the  field  of  computational  fluid  dynamics. 

First  Lieutenant  Croker  is  married  and  has  three  children. 
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